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Abstract 

Numerical relativity is the most promising tool for theoretically modeling the in¬ 
spiral and coalescence of neutron star and black hole binaries, which, in turn, are 
among the most promising sources of gravitational radiation for future detection 
by gravitational wave observatories. In this article we review numerical relativity 
approaches to modeling compact binaries. Starting with a brief introduction to the 
3+1 decomposition of Einstein’s equations, we discuss important components of 
numerical relativity, including the initial data problem, reformulations of Einstein’s 
equations, coordinate conditions, and strategies for locating and handling black 
holes on numerical grids. We focus on those approaches which currently seem most 
relevant for the compact binary problem. We then outline how these methods are 
used to model binary neutron stars and black holes, and review the current status 
of inspiral and coalescence simulations. 

Key words: 


Contents 


1 

Introduction 

2 

Decomposing Einstein’s Equations 


2.1 

Foliations of Spacetime 


2.2 

The ADM equation^ 

2.3 

Electrodynamics 

2.4 

An Illustration in Spherical Symmetry 


3 Constructing Initial Data 


3 

6 

6 

10 

15 

16 
18 


Preprint submitted to Elsevier Science 


4 February 2008 















18 


3.1 Conformal Decompositions 


3.2 The Conformal Transverse-Traceless Decomposition 


3.3 The Thin-Sandwich Decomposition 


4 Rewriting; the ADM evolution equations 


4.1 Rewriting Maxwell’s equations 


4.2 Hyperbolic formulations 


4.3 The BSSN formulation 


5 Choosing Coordinates 


5.1 Geodesic Slicing 


5.2 Maximal Slicing 


5.3 Harmonic Coordinates and Variations 


5.4 Minimal Distortion and Variations 


6 Locating Black Hole Horizons 


6.1 Locating Apparent Horizons 




6.2 Locating Event Horizons 


7 Binary Black Hole Initial Data 



7.1 The Bowen-York Approach 


7.2 The Thin-Sandwich Approach 



7.3 Comparison and Discussion 



8 Dynamical Simulations of Binary Black Holes 


8.1 Singularity Avoidance and Black Hole Excision 


8.2 Evolution of Binary Black Holes 


9 Binary Neutron Star Initial Data| 


9.1 Hydrostatic Quasi-Equilibrium 


9.2 Corotational Binaries 


20 

23 

26 

27 

29 

32 

37 

38 

39 

41 

44 

48 

49 

55 

57 

58 

64 

66 

68 

68 

74 

77 

77 

81 


2 





























9.3 Irrotational Binaries 


85 


TO Dynamical Simulations of Binary Neutron Stars 


10.1 Relativistic Hydrodynamics 


10.2 The Wilson-Mathews approximation 


10.3 Fully Self-consistent Simulations 


10.4 The Quasi-Adiabatic Inspiral of Binary Neutron Stars 


11 Summary and Outlook 


Acknowledgements 


A Notation 


B Solving the Vector Laplacian 



C Conformally Flat or Not? 



91 

91 

96 

97 
103 
106 

108 

108 

109 

111 


1 Introduction 


Promising to open a new window to the universe, a new generation of laser 
interferometer gravitational wave detectors will soon search for gravitational 
radiation. The Japanese instrument TAMA is already in operation ( [Tagoslii| 
et. alj 120011) , and the construction of the two LIGO sites in the US and the 


European instruments GEO and VIRGO is well advanced. Among the most 
promising sources for these detectors are the inspiral and merger of compact 
binaries, i.e. binaries of black holes and neutron stars. Even for these sources, 
the signal strength is expected to be much less than the detectors’ noise, so 
that sophisticated data analysis techniques will be required to extract the sig¬ 
nal from the noise (e.g. |Outler et.al . | ( |1993|) ). One such technique is matched 
filtering, in which the detector output is cross-correlated with a catalog of the¬ 
oretically predicted waveforms. The cross-correlation between a signal present 
in the data and a member of the catalog allows observers to detect signals 
that are otherwise overwhelmed by noise. Clearly, the chances of detecting a 
generic astrophysical signal depend critically on the size and quality of the 
signal catalog flElanagan and Hughcs|, |1998|). The success of the new gravita¬ 


tional wave interferometers therefore depends on accurate theoretical models 
of compact binary inspiral. 


For our purposes, the entire inspiral of compact binaries can be separated into 
three different phases. By far the longest phase is the initial quasi-equilibrium 
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inspiral phase, during which the separation between the stars decreases adi- 
abatically as energy is carried away by gravitational radiation. As the sep¬ 
aration decreases, the frequency and amplitude of the emitted gravitational 
radiation increases. Since gravitational radiation tends to circularize binary 
orbits, and since we will be mostly interested in binaries at very small sepa¬ 
ration, it is reasonable to focus on quasi-circular orbits. These quasi-circular 
orbits become unstable at the innermost stable circular orbit (ISCO), where 
the inspiral gradually enters a plunge and merger phase. The merger and coa¬ 
lescence of the stars happens on a dynamical timescale. The final stage of the 
evolution is the ringdown phase, during which the merged object settles down 
to equilibrium. 


Different techniques are commonly employed to model the binary in the dif¬ 
ferent phases. The early inspiral phase, for large binary separations, can be 
modeled very accurately with post-Newtonian methods ( [Blanchet et.al . | , |1995| ; 
parnour, Jaranowski and Schafer], |2000|, and references therein). It is gen¬ 


erally accepted that the plunge and merger phase has to be simulated by 
means of a fully self-consistent numerical relativity simulation. During the 
late ringdown phase, the merged object can be approximated as a distorted 
equilibrium object, so that perturbative techniques can be applied (Price and| 
Pullin| , |1994| : [Baker et.al\ [2001| , |2002|) . In addition to simulating the dynamical 
plunge and merger phase, numerical relativity may also be required for the late 
quasi-adiabatic inspiral phase, just outside of the ISCO, where finite-size and 
relativistic effects may become large enough for post-Newtonian point-mass 
techniques to break down. 


A number of review articles have recently appeared on various aspects of nu¬ 
merical relativity, including [D'lmci] (|2001|) on numerical relativity in general, 
Cook| ( |2000| ) on initial data, [Marti and Mullcr| (|1999Q and |Font| ( |2000| ) on nu¬ 
merical hydrodynamics in special and general relativity, [Rasio and Shapi~ro| 
( |1999|) on the coalescence of binary neutron stars, | Rcula| (|1998|) on hyper¬ 
bolic methods of Einstein’s equations, and New| ( |2002| ) on the generation of 
gravitational waves from gravitational collapse. While we will refer to these 
articles in many places, the perspective and focus of this article is very differ¬ 
ent. It is intended as a review of different numerical relativity approaches to 
the compact binary problem, and we will focus on those methods of numerical 
relativity that currently seem most promising for these purposes. Specifically, 


we will only discuss the so-called Cauchy approach in “3+T : 


i.e. 


three spa¬ 


tial dimensions (“3D”) plus timcQ, and we limit matter sources, if present, to 
perfect fluids. This means that many other promising techniques and impor¬ 
tant results of numerical relativity will not be covered, including characteristic 


1 We will provide examples of spherically symmetric configurations to illustrate 
concepts and techniques. These examples, which can be treated in 1+1 dimensions, 
are clearly marked and set in italics. 
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methods^], collisionless matter and scalar wave sources, critical phenomena, 
and results in spherical or axisymmetry. Simultaneously, we will review nu¬ 
merical relativity approaches to the compact binary problem only, and will 
not discuss post-Newtonian or perturbative methods. 


Even with this fairly focused scope of our article we will omit some important 
aspects. Perhaps most importantly, we will not cover wave extraction and 
boundary conditions for outer boundaries, since in most current applications 
in 3+1 these are only implemented fairly crudely. Some relevant references 
to more rigorous treatments include Bishop et.al .| (|1996|); [Abrahams et.al 


(|1998|) ; |Fricdrich and Nagy| ( |1999|) ; |5zilagyi et.al . | ( |2000|) ; |Calabrcsc, Lclincr 


and Tiglio| (|2002| ); [Szilagyi, Schmidt and Winicour| (|2002|) . While we will cast 


equations into a form that is suitable for numerical implementation, we will 
hardly discuss numerical methods that can be used for their solution. We 
will similarly ignore related computer science issues, including memory and 
CPU requirements and computational resources on current parallel computers. 
Some discussion of these aspects and references can be found in |Lchncr| (|2001|) . 


Loosely speaking, this article is organized into two parts. The first part, Sec¬ 
tions 2 through 6, introduces concepts and techniques of numerical relativity, 
both traditional approaches and more recent developments. The second part, 
Sections 7 through 10, reviews their applications to binary black holes and 
neutron stars. 


The first part starts by introducing the 3 + 1 decomposition and the so-called 
ADM equations in Section 2. We will see that, like Maxwell’s equations, these 
equations separate into constraint and evolution equations. In Section 3 we 
discuss strategies for solving the constraint equations and the construction of 
initial data. In Section 4 we will revisit the ADM equations, and will show that 
the evolution equations can be brought into a form that is better suited for 
numerical integrations. Before these evolution equations can be integrated, 
certain coordinate conditions have to be imposed, which we will discuss in 
Section 5. In Section 6 we review the definition of black hole horizons, and 
show how these horizons can be located in numerically generated spacetimes. 


In the second part we discuss the construction of binary black holes and neu¬ 
tron stars. Specifically, we review the initial data problem and evolution sim¬ 
ulations of binary black holes in Sections 7 and 8 and binary neutron stars in 
Sections 9 and 10. We briefly summarize in Section 11. 


Also included are three short appendixes. Appendix A explains our notation, 
Appendix B shows how the flat vector Laplace operator, which is encountered 
in several places in this article, can be solved, and Appendix C discusses some 


2 Which in fact may be a very promising alternative to the Cauchy approach even 
for the modeling of compact binaries; see, e.g. iGornez et.al . | ( |199§| ). 
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arguments which have been made in favor or against the approximation of 
spatial conformal flatness. 

As a last remark we would like to apologize to all those authors and colleagues 
whose publications we may have missed. While we have put considerable ef¬ 
fort into being as complete as possible, without doubt we have missed some 
important references. We again apologize for these unintended omissions. 


2 Decomposing Einstein’s Equations 


In this Section we briefly summarize the Arnowitt-Deser-Misner or “ADM” 
decomposition of Einstein’s equations. We will state only the most important 
relations and results and refer to the literature for more complete and rigorous 
treatments. In addition to the original article by |A rn o w itt, Deser and Misner 
(|1962|) such presentations can be found, for example, in the article by [York 
( 1979| ), the dissertations by |Smarr| (|1975|) , [fivausj (|1984|) and |Uook| ( |1990|) , and 
in 


the lecture notes of [IJaunigarte, Shapiro and Abrahams| (|1998|). 


2.1 Foliations of Spacetime 


The unification of space and time into spacetime is central to general rela¬ 
tivity and is one of its greatest aesthetic appeals. For a numerical treatment, 
however, and similarly for many mathematical treatments, it is more desirable 
to reverse this unification and recast general relativity into a so-called “3+1” 
formulation, in which a time coordinate is explicitly split from three spatial 
coordinates. Putting it differently, the four-dimensional spacetime is “carved 
up” into a family of three-dimensional spatial “slices”. 

In more technical terms, we assume that the spacetime (M, g a b) can be foliated 
into a family of non-intersecting, spacelike, three-dimensional hypersurfaces £. 
At least locally, these timeslices £ form level surfaces of a scalar function t, 
which we will later identify with the coordinate time. The 1-form 

n = dt (2.1) 


is obviously closed (df2 = 0) and has the norm 


VI 2 = = -a- 2 , 


( 2 . 2 ) 


where the lapse function a is strictly positive (see Appendix A for a summary 
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of our notation). This implies that the surfaces E are spacelike. We now define 
the timelike unit normal vector n a as 


n" = - —ag all Vbt. 


(2.3) 


Here the negative sign has been chosen so that n a points into the direction of 
increasing t. The four-dinrensional metric g ab now induces the spatial metric 

lab = 9ab + n a n b (2.4) 


on the hypersurfaces E. 

Any four-dimensional tensor can now be decomposed into spatial parts, which 
live in the hypersurfaces E, and timelike parts, which are normal to E and 
hence aligned with n a . The spatial part can be found by contracting with the 
projection operator 

7 a b = 9 ac lcb = g a b + n a n b = S a b + n a n b , , (2.5) 

and the timelike part by contracting with 

N a b = -n a n b . (2.6) 


We now define the three-dimensional covariant derivative of a spatial tensor 
by projecting all indices of a four-dimensional covariant derivative into E; for 
example 

D„T\ = (2.7) 


It is easy to show that this derivative is compatible with the spatial metric, 
Dalbc = 0, as it is supposed to be. The three-dimensional covariant derivative 
can be expressed in terms of three-dimensional connection coefficients, which, 
in a coordinate basis, are given by 

Tfoc = 2 ' l ad {ldb,c + 7 dc,b - 76c,d)- (2-8) 


The three-dimensional Riemann tensor associated with 7 ^ is defined by re¬ 
quiring that 

2 D [a D b] w c = R d cba w d R d cba n d = 0 (2.9) 
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for every spatial Wd- In terms of coordinate components, Riemann can be 
computed from 


p d 
rt abc 


prf _ p^ I pe pd _ pe pd 

L ac,b L be,a L ac L eb L bc L ea • 


( 2 . 10 ) 


Contracting the Riemann tensor yields the three-dimensional Ricci tensor 
R ab = R c acb and the three-dimensional scalar curvature R = R a a . 

It is intuitive that casting Einstein’s equations into a 3 + 1 form will neces¬ 
sitate expressing the four-dimensional Riemann tensor ^ R d abc in terms its 
three-dimensional cousin R d abc . It is also clear, however, that the latter can¬ 
not contain all the relevant information: R d abc is a purely spatial object (and 
can be computed from spatial derivatives of the spatial metric alone), while 
R d abc is a spacetime creature which also contains time-derivatives of the 
four-dimensional metric. This missing information is expressed by a tensor 
called extrinsic curvature , which describes how the slice £ is embedded in the 
spacetime M. The extrinsic curvature can be defined as 

I<ab = -7a C 76 d V (c n d) . (2.11) 


Note that K ab is spatial and symmetric by construction. An alternative ex¬ 
pression can be found in terms of the “acceleration” of normal observers 
a a = n b V b n a , 


Rab Ka^b 'bba^ J b' 


( 2 . 12 ) 


Since a a n a = 0, we immediately find for the trace of the extrinsic curvature 


K = g ab K ab = -W a n a . 


(2.13) 


In yet another equivalent expression, the extrinsic curvature can be written in 
terms of the Lie derivative of the spatial metric along the normal vector n a 


Kab = 



(2.14) 


The Lie derivative in the above equation may be thought of as the geomet¬ 
ric generalization of the partial time derivative d t . Introductions to the Lie 
derivative can be found, for example, in |Selmt \ ( |1980|) ; |Wald| ( |1984j) ; P' Inver no 
); [Baumgartc, Shapiro and Abrahams| (|1998|). Formally, the Lie deriva¬ 


tive along a vector held X a measures by how much the changes in a tensor held 
along X a differ from a mere inhnitesimal coordinate transformation generated 























by X a . For a scalar /, the Lie derivative reduces to the partial derivative 


£ x / = X b D b f = X b d b f- (2.15) 

for a vector held v a the Lie derivative is the commutator 

C x v a = X b D b v a - v b D b X a = [X, v] a , (2.16) 


and for a 1-form cu a the Lie derivative is given by 
= X b D b u a + u b D a X b . 


(2.17) 


It then follows that for a tensor T a b of rank Q) the Lie derivative is 

Cx_T a h = X c d c T a b - T\d c X a + T a c d b X c . (2.18) 

Generalization to tensors of arbitrary rank follows naturally. It is easy to verify 
that in all of the above expressions for the Lie derivative one may replace the 
partial derivative with a covariant derivative, since all connection coefficients 
cancel each other. 

Equation (2.14) most clearly illustrates the interpretation of K ab as a time- 
derivative of the spatial metric. It is therefore not surprising that spatial pro¬ 
jections of ^ R a bcd will involve the extrinsic curvature and its time derivative. 

Given its symmetry properties, ^R a bcd can be projected in three different 
ways. Projecting all four indices into E yields, after some manipulations, 
Gauss’ equation , 

Rabcd + K ac K bd - K ad K bc = Y a l q b YcYd { 4 ) R Pq rs- (2.19) 

Three spatial projections and a contraction with n a yields the Codazzi equation 
D a K bc - D b K ac = l r b l P a l q c n siA) Rr Pq s- (2.20) 


Finally, two spatial projections and two contractions with n a yield Ricci’s 
equation 

C n K ab = nV 7 9 a7 V 4) iW - ~D a D b a - K\K ac . (2.21) 

a 

In this equation the derivative of the lapse entered through the identity 

a a = D a In a. ( 2 . 22 ) 
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2.2 The ADM equations 


In the last Section we simply recorded geometrical identities relating the ge¬ 
ometry of the three-dimensional hypersurfaces E to the geometry of the em¬ 
bedding four-dimensional spacetime M. According to general relativity, the 
geometry of the latter is governed dynamically by Einstein’s equation 

^Gab = {A) Rab - \ {A) Rg ab = 8vr T ab , (2.23) 

where ^ G ab is the Einstein tensor and T ab the stress-energy tensor. We will 
now take projections of Einstein’s equations into E and n a and will use the 
Gauss, Codazzi and Ricci equations to eliminate the four-dimensional Ricci 
tensor ^ R a b■ The result will be the ADM equations, which relate three- 
dimensional curvature quantities to projections of the stress-energy tensor. 
One relation that is very useful in these derivations is 

D a V b = la c V c V b + K ac V c n b , (2.24) 


which holds for any spatial vector V a . 

Contracting Gauss’ equation (2.19) twice with the spatial metric and insert¬ 
ing (2.4) yields 

2 n a n b (4) G a6 = R + K 2 - K ab K ab . (2.25) 


We now define the total energy density as measured by a normal observer n a 
as 


p = n a n b T ab 


(2.26) 


and, using Einstein’s equation (2.23), End the Hamiltonian constraint 


R + K 2 - K ab K ab = 16tt p. 


(2.27) 


We can similarly contract the Codazzi equation (2.20) once to find the mo- 
menum constraint 


D b K b a - D a K = 8nj a , 


(2.28) 


where 

ja = -rf a n c T bc 


(2.29) 
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is the momentum density (mass current) as measured by a normal observer n a . 
The Hamiltonian and momenum constraints are called constraint equations 
because they only involve spatial quantities and their spatial derivatives. They 
therefore have to hold on each individual spatial slice E - in fact they are 
the necessary and sufficient integrability conditions for the embedding of the 
spatial slices (E, 7 a b,K a b) in the spacetime (. M,g a b ). 

Evolution equations that describe how data y a b and K ab evolve in time, from 
one spatial slice to the next, can be found from equation (2.14) and the Ricci 
equation (2.21). However, the Lie derivative along n a , £ n , is not a natural time 
derivative orthogonal to the spatial slices, since n a is not dual to the surface 
1 -form i.e. their dot product is not unity but rather 

n a Vt a = —ag^Vb^V Cl t = a -1 . (2.30) 


Instead, the vector 
t a = an a + (3 a 


(2.31) 


is dual to f2 for any spatial shift vector j3 a . The lapse a and the shift /3 a together 
determine how the coordinates evolve from one slice E to the next. The lapse 
determines how much proper time elapses between timeslices along the normal 
vector n a , while the shift determines by how much spatial coordinates are 
shifted with respect to the normal vector. 

The Lie derivative along t a is a natural time derivative, because the duality 
t a Q a = t a V a t = 1 (2.32) 


implies that the integral curves of t a are naturally parametrized by t. As a 
consequence, all (infinitesimal) vectors t a originating on one spatial slice Ei 
will end on the same spatial slice E 2 (unlike the corresponding vectors n a , 
which generally would end on different slices). This also implies that the Lie 
derivative of any spatial tensor along t a is again spatial (see, e.g., Problem 
8.14 in |Lightman et.al . | (|1975|) for an illustration). 


Rewriting equation (2.14) in terms of t a yields the evolution equation for the 
spatial metric 


'iab 2 OlK a b T 


(2.33) 
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The evolution equation for the extrinsic curvature can be found by combining 
Ricci’s equation (2.21) with Einstein’s equations (2.23) 


£t kab — ~D a D b a + a(R ab — 2 K ac K c b + RK ab ) 
a8n(S a b 274 6 ( p)) T RpRabi 


(2.34) 


where is the spatial projection of the stress-energy tensor 


Safe = laclMT Cd , 


(2.35) 


and S' its trace S = 7 a 6 S a 6 . 

While the two constraint equations (2.27) and (2.28) constrain 7 afe and /d a6 
on every spatial slice E, the evolution equations (2.33) and (2.34) describe 
how these quantities evolve from one slice to the next. It can be shown that 
the evolution equations preserve the constraint equations, meaning that if the 
constraints hold on one slice, they will continue to hold on later slices. This 
structure is very similar to that of Maxwell’s equations, as we will discuss in 
Section 2.3. 

So far, we have made no assumptions about the choice of coordinates, and 
have expressed all quantities in a coordinate-independent way. It is quite in¬ 
tuitive, though, that things will simplify if we adopt a coordinate system that 
reflects our 3+1 split of spacetime. To do so, we introduce a basis of spatial 
vectors e, (where i — 1 , 2 , 3, see Appendix A for a summary of our notation) 
that span each slice E, so that f2 a (ej) a = 0. It can be shown that this condition 
is preserved if the spatial vectors are Lie dragged along t a . As the fourth basis 
vector we pick (eo) a = £ a , or, in components, t a = ( 1 , 0 , 0 , 0 ). As an immediate 
consequence, the Lie derivative along t a reduces to a partial derivative with 
respect to t. Since n l = n a (ej) a = af2 a (ej) a = 0, the covariant, spatial com¬ 
ponents of the normal vector vanish. Since contractions of any spatial tensor 
with the normal vector are zero, this implies that the zeroth components of 
contravariant spatial tensors have to vanish. For the shift vector, for example, 
/3 a n a = P°n 0 = 0, so that we can write 

(3 a = (0,/T). (2.36) 


Solving (2.31) for n a then yields the contravariant components of the shift 
vector 


n“ = 1 ( 1 , -fi) 
a 


(2.37) 
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(2.38) 


and, since n a n a — — 1 , 
n a = (-a, 0 , 0 , 0 ). 


From the definition of the spatial metric (2.4) we find 7 ^ = g^. Since zeroth 
components of contravariant, spatial tensors are zero, we also have 7 “° = 0 . 
The inverse metric can therefore be expressed as 


g ab = 


( -a .- 2 


cC 2 fi l 

\a~ 2 p j 7 ij -a~ 2 f3 i f3 j 


(2.39) 


Using n t = 0, 7 ^ = g l3 and 7 “° = 0 it is also possible to show that 7 ^ and 
7 y are inverses r f kr )kj — j, so that they can be used to raise and lower 
indices of spatial tensors. Inverting (2.39) we find the components of the four¬ 
dimensional metric 


9ab 


/ -a 2 + /3 k /3 k 
\ Pj 



(2.40) 


and equivalently the line element 

ds 2 = —a 2 dt 2 + 7 ij(dx l + (3 l dt)(dx J + ftdt). (2-41) 


The entire content of a spatial tensor is encoded in its spatial components 
alone. We can therefore restrict the constraint and evolution equations to 
spatial components. Moreover, since in our coordinate system the zeroth com¬ 
ponents of contravariant, spatial tensors are zero, we can also restrict all con¬ 
tractions to spatial components. The connection coefficients (2.8) reduce to 

= 2 ^( 7 ij,k + hk ,3 - 7 jk,i), (2.42) 


and expressing the Ricci tensor (2.10) in terms of second derivatives of the 
metric yields 


R ij 2 


7 ( 7 kj,il "b 'Yil,kj ''ikl.i/j 'Yij,klj 

+ 7 fc/ (r™r mfci - r™r, mkl 


(2.43) 


With these simplifications, the Hamiltonian constraint (2.27) now becomes 


R + K 2 - KijK ij = 16vrp, 


(2.44) 
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the momenum constraint 


DjK\ - DjK - 87 rji, 


(2.45) 


the evolution equation for the spatial metric (2.33) 


dtlij — —ZatKij + D{f3j + Djfii, 


(2.46) 


and the evolution equation for the extrinsic curvature (2.34) 


dtKjj — — D,DjQ + a(Rij — 2 Ki k K k ^ + KKij ) 

-a8n(Sij - ij(S - p)) (2.47) 

+P k D k K ij + K ik Dj(3 k + K kj Di(3 k . 


The shift terms in the last two equations arise from the Lie derivatives £^ 7 ^ 
and CpKij. In (2.46) it is convenient to express the Lie derivative in terms of 
the covariant derivative Di which eliminates the term /3 k D k ^ij, but in (2.47) 
the covariant derivative may be replaced with the partial derivative di in these 
terms. 


Equations (2.44) to (2.47) are commonly referred to as the ADM equations^]. 
As it turns out, these equations are not yet in a form that is generally well 
suited for numerical implementation. In Section 4 we will see that the stabil¬ 
ity properties of numerical implementations can be improved by introducing 
certain new auxiliary functions and rewriting the ADM equations in terms 
of these functions. Most current numerical relativity codes are based on such 
reformulations of the ADM equation. 

Note that the ADM equations only determine the spatial metric 7 ^ and the 
extrinsic curvature K,j , but not the lapse a or the shift (3 l . The latter deter¬ 
mine how the coordinates evolve from one timeslice to the next and reflect 
the coordinate freedom of general relativity. Choosing coordinates that are 
suitable for the situation that one wishes to simulate is central to the success 
of the simulation. I 11 Section 5 we will discuss strategies for both choosing and 
numerically implementing various coordinate conditions. 


For later purposes it is useful to take the traces of the evolution equations 
(2.46) and (2.47). Since c^lny = Y j dt”/ij, we find 


d t I 117 1 / 2 = —aK + Di/3 l , 


(2.48) 


3 Arnowitt, Deser and Misneij Q1962 ) originally derived these equations in terms of 
the conjugate momenta nij instead of the extrinsic curvature. 
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for the determinant of the metric 7, and combining the Hamiltonian constraint 
(2.44) with the trace of (2.47) yields 

d t K = —D 2 a + a (k^K ij + 4vr (p + S)) + fiDiK. (2.49) 

Here D 2 = 7 U AA is the covariant Laplace operator associated with 7^. 


2.3 Electrodynamics 


To discuss the structure of the ADM equations it is very instructive to compare 
them to the equations of electrodynamics. Maxwell’s equations in flat space 
naturally split into a set of constraint equations 


D i E i = 4np e (2.50) 

D i B i = 0 (2.51) 

which hold at each instant of time, and a set of evolution equations 

d t Ei = e ijk D>B k - 47 iJi (2.52) 

d t B l = -e ijk D^E k (2.53) 


which describe the evolution of the electric held A and the magnetic held A 
from one instant of time to the next. Here p e and J* are the charge density and 
current. Note that the evolution equations preserve the constraints so that, if 
they are satished at any time, they are automatically satished at all times. 

It is often useful to introduce the vector potential A a = (—<3>, Ai) and write 
B l as a curl of Ai 

B i = e ijk DjA k (2.54) 

so that B l satishes the constraint (2.51) automatically. Maxwell’s equations 
can now be rewritten as two evolution equations for A* and E t 


d t Ai = —E { - AT (2.55) 

dtEi = —D^DjAi + D t D 3 A j — 47rJj (2.56) 


together with the constraint equation (2.50) (see, e.g., |Jackson| (|1975| )). The 


gauge quantity T is, like the lapse and shift in the ADM equations, undeter¬ 
mined by the equations, and has to be chosen independently. 
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Interestingly, the evolution equations (2.55) and (2.56) are quite similar to the 
ADM evolution equations (2.46) and (2.47), which can be seen by identifying 
Ai with jij and E, with K io . Both right hand sides of (2.46) and (2.55) involve 
a held variable and a derivative of a gauge variable, while both right hand 
sides of (2.47) and (2.56) involve matter sources as well as second derivatives 
of the second held variable (which in (2.47) are hidden in the Ricci tensor 
(2.43)). In Section 4.1 we will further explore these similarities. 


2-4 An Illustration in Spherical Symmetry 


To illustrate the concepts introduced in this Section, it is useful to work out 
some of the expressions in spherical symmetry. Throughout this Article we 
will return to this example. 

Example 2.1 The general form of a spherically symmetric spacetime metric 
is 

ds 2 = —Adt 2 + 2 Bdtdr + Cdr 2 + D(dd 2 + sin 2 Odf 2 ), (2.57) 


where the coefficients A, B, C and D are functions of time t and radius r 
only (see, e.g., equation (14-29) in \D dnverndj (\1992) ). We now introduce a 
new radial coordinate by the transformation 

r -> r = ( D/C) 11 2 , (2.58) 


which, after dropping all tildes, brings the metric into the isotropic form 


ds 2 = —Adt 2 + 2 Bdtdr + C{dr 2 + r 2 d 2 Q), 


(2.59) 


where d 2 D = dO 2 + sin 2 9d<f 2 . Comparing with the metric (2-41), we can now 
identify B with (3 r and, since /3 r has to be the only non-vanishing component of 
the shift vector, A with a 2 — /3 r /3 r . Following convention, we will also rewrite C 
as 0 4 , where -0 is called the “conformal factor” (compare Section 3). Defining 
/3 as the contravariant component (3 r , we have 

ds 2 = —(a 2 — il) A (I 2 )dt 2 + 2 ^fldrdt + -0 4 (dr 2 + r 2 d 2 fl). (2.60) 
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The only non-vanishing, spatial connection coefficients can now be found to be 


r 


r 

rr 



= 2 — 
sin 2 9 T r ee 


~p0 _ ~p0 _ ~P</ > _ I ^ 

1 re - 1 0r - 1 r<t> ~ 1 <j>r ~ Z ^ 


T% = — sin 9 cos 9 
it, = r*„ = cote. 


(2.61) 


Here and in all following examples a prime denotes partial derivative with 
respect to radius r, and a dot will denote partial derivative with respect to 
time t. As expected, the connection reduces to that of a flat metric in spherical 
symmetry for if = 1. The non-vanishing components of the extrinsic curvature 
can be computed from equation (2.f6) according to 


K r = _ 2 (t _ 

r ^ 2) 

Kd e = K \ = --[l-^ 

v a Kip ip 


(2.62) 


and the trace of the extrinsic curvature is 

aK = ~(ip - flifl) + 4 (fir 2 )'- 
ip r z 


(2.63) 


Note that subtracting the two equations in (2.62) yields an equation for the 
shift 


r/3' -13 = a(K r r - K%) r 


(2.64) 


which can be integrated to yield 


(3 = r 


a (A"; - K%) 
r 


dr. 


(2.65) 


It may seem surprising that we fijid an equation for the shift, even though 
we have emphasized on several occasions that Einstein’s equations do not de¬ 
termine the lapse and the shift, and that, representing the coordinate freedom 
in general relativity, they have to be chosen independently. However, we have 
already used up the spatial coordinate freedom in equation (2.58), by bringing 
the metric into the isotropic form (2.60). The shift condition (2.65) ensures 
that the metric remains isotropic for all times, and it can therefore no longer 
be chosen freely. The lapse, however, is still undetermined. 
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One could carry the exercise further, and derive, for example, the Ricci ten¬ 
sor Rij from the connection (2.61), which would be straight-forward but quite 
lengthy. Instead we will postpone this until we develop a formalism in Section 3 
that will simplify this exercise significantly. 


3 Constructing Initial Data 


In this Section we will present strategies for solving the constraint equa¬ 
tions (2.44) and (2.45). Since a review of the initial value problem for nu¬ 
merical relativity has appeared very recently ( |CJool4 P000| ), we will focus only 
on the the most important results that are often employed for the construction 
of binary black holes and neutron stars. In this Section we will review the most 
common formalisms, and will defer the discussion of particular solutions for 
binary black holes and neutron stars to Sections 7 and 9. Parts of this Section 
are based on the lecture notes of |Baumgartc, Shapiro and Abraham^ (|1998|). 


3.1 Conformal Decompositions 


Most approaches to solving the relativistic initial value problem involve a 
conformal decomposition, in which the physical metric 7^ is written as a 
product of a conformal factor if and an auxiliary metric, usually referred to 
as the conformally related or background metric 7^, 

7 ij = ij (3-1) 


( |Lichnerowicz| , |1944| ; |Yorl4 |1971| , |1972|) . Taking ^ to the fourth power turns 
out to be convenient, but is otherwise arbitrary. For other purposes it is advan¬ 
tageous to consider conformal transformations of the four-dimensional space- 
time metric g a b, but here we will only consider conformal transformations of 
the spatial metric 7^. 


Superficially, we can think of the decomposition (3.1) as a mathematical trick, 
namely writing one unknown as a product of two unknowns, which makes the 
solving of some equations a little easier. At a deeper level, the conformal 
transformation (3.1) defines an equivalence class of manifolds and metrics, 
which are all related by the conformal metric 7y = 7~ 1 / 3 7p-, where 7 is the 
determinant of the metric 7^. I 11 this natural definition 7 = 1, but other 
normalizations can be chosen. A metric that is conformally related to the flat 
spatial metric, 7^ is called conformally flat. 
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(3.2) 


From (3.1), the connection coefficients transform according to 
r ) k = f ) k + 2 (S i j d k In ip + 5\dj hup - 7 jk f l di hup), 


from which it is easy to verify that A 7 ij = 0. For the Ricci tensor (2.43) we 
find 

Rij = Rij - 2 (A A In ip + %*/ lm DiD m In 0 ) 

+4 ((A In 0)(A In 0) - 7^7 Zm (A In 0)(An In 0)), 


and for the curvature scalar 

R = xp~ A R - 8 ip~ b D 2 ip. 


(3.4) 


Here D 2 = 7 tJ I); D,, Rij and R are the covariant Laplace operator, Ricci tensor 
and scalar curvature associated with 7 ^. R,j can be computed by inserting 
into equation (2.43). 

Example 3.1 Returning to Example 2.1, we now see that the spatial part of 
the metric (2.60) is a conformal factor 0 4 times the flat metric rjij in spherical 
coordinates. Since any spherically symmetric metric can be brought into the 
form (2.60) without loss of generality, we have shown that any spherically 
symmetric metric is (spatially) conformally flat. 

Instead of computing curvature quantities from the spatial metric jij, it is 
now much easier to employ the above formalism. The connection coefficients 
(2.61), for example, can be found by adding the connection coefficients T l - k 
of a flat metric in spherical coordinates to spatial derivatives of 0 according 
to (3.2). Since 7 ^ = r/ y - is flat, the conformally related Ricci tensor vanishes 
Rij = 0 and the physical Ricci tensor is given by the derivatives of ip in (3.3) 
alojie. The Ricci scalar, fi7ially, reduces to 

R = -8 0- 5 L> 2 A (3.5) 


a remarkably simple expression. 

It is also convenient to split the extrinsic curvature K tJ into its trace K and 
a traceless part A t j 

Kij = Aij + ^ 7 ijK. (3.6) 

Solving the initial value problem usually proceeds by further decomposing the 
traceless part A,j, which is done slightly differently in different approaches. 
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In the following we will briefly discuss the transverse-traceless and the thin- 
sandwich decomposition. 


3.2 The Conformal Transverse-Traceless Decomposition 


In the conformal transverse-traceless decomposition, we first introduce the 
conformal traceless extrinsic curvature by defining 

A ij = fj~ 10 A ij (3.7) 


',-2 


(see, e.g., |York| (|1979|) ). Other choices for the 


and accordingly A i3 = if 
exponent of the conformal factor are possible as well - and in fact we will use 
a different scaling in Section 4.3 - but for our purposes here the exponent —10 
is particularly convenient since any symmetric traceless tensor S'*- 7 satisfies 
DjS ij = r 10 /;,((.• "'S'-'). 

Any symmetric, traceless tensor can be split into a transverse-traceless tensor, 
which is divergenceless, and a longitudinal part, which can be written as a 
symmetric, traceless gradient of a vector (|York| , |T973|) . We can therefore further 
decompose A*- 7 as 


Aij _ A l J I A l J 

IX - IX'Jf'JT' IX , 


(3,8) 


where the transverse part is divergenceless 
D 3 Af T = 0 


(3.9) 


and where the longitudinal part satisfies 

A i l = D { W j + D j W { --f j D k W k = ( LW) ij . 

3 


(3.10) 


Here W is a vector potential, and it is easy to see that the longitudinal operator 
or vector gradient L produces a symmetric, traceless tensor. The divergence 
of A lj can be written as 


DjA ij = D J A i l = D 3 (LW) ij = D 2 W i + \D\DjW j ) + If fV' 

= (A l W)*, 


(3.11) 


where is the vector Laplacian. 
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Inserting the conformally related quantities into (2.44) now yields the Hamil¬ 
tonian constraint 


8 D 2 if 


ipR - -i) b K 2 + if 


-7 


AijA V = —1677-0 p, 


(3.12) 


while the momentum constraint (2.45) becomes 

(AlWY - -if^DjK = 8vr0 10 /. (3.13) 

3 

Example 3.2 As a consequence of the Bianchi identities, not all of Einstein’s 
equations are independent. This redundancy can be exploited in numerical evo¬ 
lution calculations, where some quantities can be solved for using either con¬ 
straint or evolution equations. The conformal factor if, for example, satisfies 
an evolution equation (which in spherical symmetry can be found by solv¬ 
ing equation (2.63) for if) and the Hamiltonian constraint (3.12). In fact, 
the evolution equations can be completely eliminated in spherical symmetry, 
which reflects the fact that in spherical symmetry the gravitational fields do 
not carry any dynamical degrees of freedom. In constrained evolution codes 
some of the evolution equations are replaced by constraint equations, while in 
unconstrained evolution codes the gravitational fields are computed from the 
evolution equations alone, while the constraint equations are typically used to 
monitor the quality of the numerical solution. Constrained evolution has been 
very popular in spherical and axial symmetry. In full 3D, however, solving 
the elliptic constraint equation is very expensive computationally, so that most 
current codes in 3 + 1 use unconstrained evolution. 

Given a solution to the constraint equations, one can find the mass, linear and 
angular momenta from the asymptotic behavior of the solution. For asymp¬ 
totically flat solutions, the total (ADM) energy is 

M — “ 2 tt / D K fd 2 Si, (3.14) 

OO 

the linear momentum is 

P i = — cf K ij d 2 Sj, (3.15) 

877 J 


and, in Cartesian coordinates, the angular momentum is 
A = — x j K kl d 2 S t 

* 877 7 

OO 


(3.16) 
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(see, e.g., P Murchadha and York| (|1974j) ; [Bowen and York| ( |1980| ); |York 

«))• 


At this point it is instructive to count degrees of freedom. We started out 
with six independent variables in both the spatial metric 7 ^ and the extrinsic 
curvature K t] . Splitting off the conformal factor 0 left five degrees of freedom 
in the conformally related metric 7 ^ (specifying 7 ). Of the six independent 
variables in K lJ we moved one into its trace K, two into A l ^ T (which is sym¬ 
metric, traceless, and divergenceless), and three into A l [ (which is reflected 
in its representation by a vector). Of the twelve original degrees of freedom, 
the constraint equations determine only four, namely the conformal factor 0 
(Hamiltonian constraint) and the longitudinal part of the traceless extrinsic 
curvature (momentum constraint). Four of the remaining degrees of freedom 
are associated with the coordinate freedom - three spatial coordinates hidden 
in the spatial metric and one that determines the evolution in time that is 
often identified with K. This leaves four physical degrees of freedom undeter¬ 
mined - two in the conformally related metric 77 ,-, and two in the transverse 
part of traceless extrinsic curvature A l ^ T . These two freely specifiable degrees 
of freedom carry the dynamical degrees of freedom of the gravitational fields. 
All others are either fixed by the constraint equations or represent coordinate 
freedom. 

It is obvious that the background data 77 ,-, A l ^ T and K have to be chosen before 
the constraints can be solved for 0 and Al[ . The choice of background data has 
to be made in accordance with the physical or astrophysical situation that one 
wants to represent. Physically, the choice reflects the amount of gravitational 
wave content in the background. In many situations, as for example for close 
compact binaries, it is not a priori clear what choices correspond to the desired 
astrophysical scenario and reflect the past inspiral history. Arguments have 
then be made that the background data can be approximated reasonably well 
by conformal flatness 77 , = 77 ^, maximal slicing K = 0 (see Section 5.2) and 
Alp T = 0. While this choice is somewhat controversial (see the discussion in 
Appendix C), it has the great benefit of dramatically simplifying the constraint 
equations (3.12) and (3.13). 

It is quite remarkable that the momentum constraint (3.13) is a linear equa¬ 
tion for the vector potential W\ For maximal slicing K = 0 and in vacuum 
j — 0 it also decouples from the Hamiltonian constraint, and can be solved 
independently of a solution for 0. For conformal flatness, the operator sim¬ 
plifies significantly (see also Appendix B), and in fact we will discuss analytic 
solutions in Section 7.1. 

Assuming conformal flatness, the Hamiltonian constraint (3.12) also simplifies 
significantly, since D 2 reduces to the flat Laplace operator and R = 0. Maximal 
slicing further simplifies the equation. 
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Example 3.3 As we have seen in Example 3.1, the metric (2.60) is confor¬ 
mally flat. If we are further interested in time-symmetric vacuum solutions 
(whereby K t] = —Kij = 0 and p = 0/, the Hamiltonian constraint (3.12) 
reduces to the simple Laplace equation 

A fla V = 0 (3.17) 


in spherical symmetry (where A flat is the flat Laplace operator). Assuming 
asymptotic flatness, so that —> 1 as r —> oo, we can write the solution as 




m 

1 + Tr 


(3.18) 


where m is a constant. Evaluating the ADM mass (3.1 f) of this solution shows 
that M — m. We have therefore rediscovered the spatial part of the Schwarz- 
schild solution in isotropic coordinates. We will see in Example 6.1 that the 
surface r = m /2 is the event (and apparent) horizon of a black hole. 


It can be shown that the coordinate transformation x l —> (m/2) 2 x l /r 2 maps 
every point inside r = m/2 into a point outside r = m/2. Moreover, the 
transformation maps the metric into itself, making it an isometry. Applying 
the coordinate transformation twice yields the identity transformation, and 
points on the horizon r = m /2 are mapped into themselves. We can therefore 
think of the transformation as a reflection across r = m/2. The geometry 
close to the origin is therefore identical to the geometry near infinity, which 
suggests that we can think of the geometry as two asymptotically flat sheets, or 
else as two separate but identical universes, which are connected by a throat or 
Einstein-Rosen bridge ( \Einstem and Roser\ , \1935) . See Figure 31.5 in \Misner, 
Thorne and Wheelei\ \1973f for an embedding diagram. 


It is remarkable that equation (3.17) is linear. As an immediate consequence, 
time-symmetric initial data containing multiple black holes can be constructed 
by simply adding several terms of the form (3.18). We will return to the con¬ 
struction of binary black hole initial data in Section 7. 


3.3 The Thin-Sandwich Decomposition 


Solving the conformal transverse-traceless decomposition yields data 7 ^ and 
K i:i intrinsic to one spatial slice £, but this solution does not tell us anything 
about how it will evolve in time away from £, nor does the formalism allow 
us to determine any such time-evolution. In some circumstances, for exam¬ 
ple when we are interested in constructing equilibrium or quasi-equilibrium 
solutions, we would like to construct data such that they do have a certain 
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time-evolution. The thin-sandwich approach offers an alternative that does 
allow us to determine the evolution of the spatial metric. Instead of providing 
data for and K l3 on one timeslice, it provides data for 7 ^ on two timeslices, 
or, in the limit of infinitesimal separation of the two slices, data for 7 i3 and 
its time derivative. The original thin-sandwich conjecture goes back to |Baier- 
lcin, Sharp and Wheeler| ( |1962| ) (see also the discussion in |Misner, Thorne and 


W heeler] (|1973|)), but here we will focus on a more recent, conformal formu¬ 


lation ( |York| , |1999| ; |Cook| , |2000|) . A similar formulation has been developed 
independently (and earlier) by |Wilson and Mathcws| ( |1989| , |1995|) (see Section 
10 . 2 ). 


Following |Cook| (|2000|) , to whom the reader is referred for a more detailed 
treatment, we start by defining u^ as the traceless part of the time derivative 
of the spatial metric, 

Uij = 7 1/3 ^(7 _1/3 7ii), (3T9) 


in terms of which the evolution equation (2.46) becomes 
u ij = -2 aA ij + {L(3) ij . 


(3.20) 


Here L is the vector gradient defined in (3.10), except in terms of the physical 
metric 7 ^. Using[] 

= -0 4 u i:j (3.21) 

together with the conformal scaling (3.7) and the identity (L /?)*- 7 = ^^(Z/?)*- 7 , 
we can rewrite equation (3.20) as 

A ij = ^ ([Lj3) ij - u ij ) . (3.22) 

This equation relates A *- 7 to the shift vector /3 l . Inserting this equation into 
the momentum constraint (2.45), we can therefore derive an equation for the 
shift 


(A l ( 3) 1 — (L(3) l iD 3 liAaijj 6 ) = 

ai/j~ (i Dj(a~ 1 , ijj e u l - i) -\— aD l K + 167ra-0 4 J*- 

3 


(3.23) 


A solution of the thin-sandwich formulation can now be constructed as follows. 
We choose the background metric 7 ^ as well as its time derivative Uij. Given 

1 This follows from the identifications Uij = dt'yij and 7 l -*Uij = 0. 
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choices for the lapse a and the trace of the extrinsic curvature K , we can then 
solve the Hamiltonian constraint (3.12) and the momentum constraint (3.23) 
for the conformal factor 0 and the shift 0Q With these solutions, we can 
then construct A from (3.22) and finally the physical quantities 7 y and K\j. 

It is again instructive to count the degrees of freedom, and to compare with 
the transverse-traceless decomposition of Section 3.2. There, we found that of 
the twelve independent variables in 7 ^ and Ky, four were determined by the 
constraint equations, four were related to the coordinate freedom, and four 
represented the dynamical degrees of freedom of general relativity. The latter 
eight can be chosen freely. In the thin-sandwich formalism, we count a total 
of sixteen independent variables, of which we can freely choose twelve: five 
each in 7 ^ and Uij, and one each for a and K. The four remaining variables, 
0 and f3\ are then determined by the constraint equations. The four new 
independent variables are accounted for by the lapse a and the shift /3 l , which 
are absent from the transverse-traceless decomposition. That approach only 
deals with quantities intrinsic to one spatial slice E, and hence only requires 
coordinates on E. The thin-sandwich approach, on the other hand, also takes 
into account the evolution of the metric away from that slice, and therefore 
requires coordinates in a neighborhood of E. As a consequence, the lapse a 
and the shift /3 l , which describe the evolution of the coordinates away from E, 
appear in the thin-sandwich approach, but not in the the transverse-traceless 
decomposition. The four new free degrees of freedom hence reflect the time 
derivatives of the coordinates. 

The thin-sandwich approach is particularly useful for the construction of equi¬ 
librium or quasi-equilibrium data, for which it is natural to choose 

= 0. (3.24) 


For equilibrium data it is also natural to choose K = 0 and d t K = 0, which 
from equation (2.49) yields the maximal slicing condition 

D 2 a = a ( AijA ij + 4t r(p + S')) . (3.25) 


We will discuss maximal slicing in more detail in Section 5.2. Expressing D 2 
in terms of D 2 and using the conformal transformations of Section 3.1, we find 
that this equation can be combined with the Hamiltonian constraint to yield 


A L (a0) = cc0 -0 A^A 13 + -R + 27T0 (p + 2 S) 


(3.26) 


5 Strictly speaking, the “densitized” lapse a = 7 -1 / 2 a (see also equation (4.16) 
below) should be fixed instead of the lapse a ( [York , 199S| ). This distinction is un¬ 
necessary if the lapse is determined through maximal slicing. 
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Equations (3.12), (3.23) and (3.26) now determine the solutions for ip, j3 l and 
a. The equations further simply if we assume conformal flatness, in which case 
they reduce to 


A fla V = -\ip~ 7 A i:j A ij - 2t uP 5 p (3.27) 

8 

(A = 2A ij Dj(aip~ 6 ) + 16tt aip 4 f (3.28) 

ip~ s AijA ij + 2mp A (p + 2S)j . (3.29) 



Here A flat and A® at are the flat Laplacian and vector Laplacian. Strategies for 
solving the flat vector Laplacian will be discussed in Appendix B. Interestingly, 
we will re-discover the shift condition (3.28) in Section 5.4 and will find that 
it is identical to minimal distortion. The thin-sandwich formalism therefore 
reduces to the Hamiltonian constraint for the conformal factor, the minimal 
distortion condition for the shift, and the maximal slicing condition for the 
lapse. 

If initial data for a time evolution calculation are constructed from the trans- 
verse-traceless decomposition, then the lapse and shift have to be chosen inde¬ 
pendently of the construction of initial data. The thin-sandwich formalism, on 
the other hand, provides a lapse and a shift together with the initial data 7 ^ 
and Kij. Obviously, once the initial data are determined, the lapse and shift 
can always be chosen freely in performing subsequent evolution calculations. 
However, the original relation between the time derivative of 7 ^ and Uij only 
applies when the lapse and shift of the thin-sandwhich solution are employed. 


4 Rewriting the ADM evolution equations 


The ADM equations as presented in Section 2.2 form the basis for most 3+1 
decompositions of the Einstein equations and their numerical implementa¬ 
tions. As it turns out, the evolution equations are not yet in their most desir¬ 
able form, and straight-forward implementations in three spatial dimensions 
typically develop instabilities very quickly. 

In this Section we will discuss the draw-backs of the ADM equations and 
will present possible alternatives. We will Erst illustrate some of the relevant 
issues with an electrodynamics analogy, and will then present two different 
re-formulations of the ADM equations. 
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4-1 Rewriting Maxwell’s equations 


In Section 2.3 we showed that Maxwell’s equations, when written as the two 
evolution equations (2.55) and (2.56) 


A = ~E t - Ad) (4.1) 

d t Ei = -DWjAi + DiD j Aj - AirR (4.2) 

and the constraint equation (2.50) 

D% = 4vr p e , (4.3) 


share some of the structure of the ADM equations. In fact they also share 
some of their draw-backs. To illustrate these, take a time derivative of (4.1) 
and insert (4.2) to form a single equation for the vector potential A 

- d 2 t Ai + DWjAi ~ DW j Aj = DA $ - AnR. (4.4) 


It is obvious that this equation would be a simple wave equation for the com¬ 
ponents A if it weren’t for the mixed derivative term DiD^Aj. In general 
relativity the situation is very similar, since the Ricci tensor Rij (2.43) on the 
right hand side of equation (2.47) also contains mixed derivative terms in addi¬ 
tion to a Laplace operator acting on r ) l] . Without these mixed derivatives, the 
ADM equations could be written as wave equations for the components of the 
spatial metric and would be manifestly hyperbolic (see, e.g., |Fricdricb] Ql996|) 
for a discussion). This is unfortunate, since it would be desirable in many 
ways to deal with a hyperbolic system of equations. Mathematical theorems 
would guarantee the existence and uniqueness of solutions, and for numerical 
purposes one could bring the equations into a form that allows for the ap¬ 
plication of flux-conservative schemes that have been developed and tested in 
other fields of computational physics, and quite in general one might feel more 
comfortable that numerical implementations will produce stable evolutions. 


These considerations suggest that it might be desirable to eliminate the mixed 
derivative terms. In electrodynamics, three different approaches can be taken 
to eliminate the DjZR A term: one can make a special gauge choice, one can 
bring Maxwell’s equations into an explicitly hyperbolic form, or one can intro¬ 
duce an auxiliary variable. For the remainder of this section, we will briefly 
discuss each one of these three strategies. 


27 







The most straightforward approach is to choose a gauge such that the term 
DjD^Aj disappears. This can be achieved, for example, in the Lorentz gauge 

d t $ = —D i A i (4.5) 


for which these two terms cancel in equation (4.4), reducing it to a wave 
equation for the vector potential A t (see |Jackson| (|1975|) ). Another possibility 
is the Coulomb gauge D L Ai = 0, which results in an elliptic equation for $. 

In general relativity, an analogous approach can be taken by choosing har¬ 
monic coordinates (see Section 5.3), which bring the equations into a mani¬ 
festly hyperbolic form. This was first realized by De Ponder (1921) and Cano 


zos| ( |1922|) , and many recent hyperbolic formulations of Einstein’s equations 


are based on this gauge choice (e.g. |Choquct-Bruhat| ( |1952| , [TT)(I2| ); [Fischer and 
Marsdenj (|1972|) ). For the purpose of numerical simulations this strategy does 
not seem very promising, since harmonic coordinates may not be the opti¬ 
mal coordinate choice for the astrophysical situation at hand (but see [Landry 
and Tcukolskyl (|1999|) , where the traditional ADM approach together with 
harmonic coordinates and finely tuned finite differencing has been used to 
simulate the merger of binary neutron stars; see also |Garfinklc| (|20Q2| ) ). More¬ 
over, harmonic coordinates may develop pathologies which could prematurely 
end a numerical simulation QAlcubicrrc , 1997 ; Alcubicrrc and Masso| , 1998 ). It 
is therefore more desirable to preserve the coordinate freedom, and to adopt 
a different approach to re-writing the equations. 

An alternative, gauge-covariant approach to bringing Maxwell’s equations into 
an explicitly hyperbolic form is to take a time derivative of (4.2) instead of (4.1), 
which yields ( Abrahams and York , 1997 ) 


c^E, = DiD^-E, - Dfi) - DjD^-Et - D t i) - d t J„ 


(4.6) 


Using the constraint (4.3) we can eliminate the first term and find a wave 
equation for Ei 

- d 2 E, + DjD j Ei = d t J t + 4tt D iPe . (4.7) 


Interestingly, the gauge dependent quantities A, and $ have disappeared from 
this equation, and, quoting [Abrahams and York| ( |1997|) , “the dynamics of 
electromagnetism have been cleanly separated from the gauge-dependent evo¬ 
lution of the vector and scalar potentials.” We will discuss similar approaches 
in general relativity in Section 4.2. 


While equation (4.7) is aesthetically very appealing, it also reveals some po¬ 
tential disadvantages for simulations that involve matter. In evolution calcu- 































































lations of neutron stars, for example, the non-smoothness of the matter and 
fields on the stellar surface or across shocks always pose numerical difficulties. 
It is to be expected that these would only get worse if further derivatives of 
the matter variables have to be taken, as in equation (4.7). 

This suggests a third approach to re-writing Maxwell’s equations, namely by 
introducing an auxiliary variable 

T = D%. (4.8) 


Inserting this into (4.2) yields 

d t Ei = -PjD'Ai + DiT - AnJi. 


(4.9) 


We now elevate T’s status to that of a new independent variable, and derive 
its evolution equation from (4.1) 


dtT = d t D i A i = D%Ai = -D^i - D i D i A 
= -DiD l <& - Anp e . 


(4.10) 


Note that we have used the constraint equation (4.3) in the last equality, 
similar to how we used the same constraint to arrive at the wave equation (4.7). 
We will see in equation (4.31) below that this step is crucial for stabilizing the 
system. 

Equations (4.1), (4.9) and (4.10) are now the evolution equations in this new 
formulation, and equations (4.3) and (4.8) are the constraint equations. In this 
formulation the mixed derivative term D t D J A 3 has been eliminated without 
using up the gauge freedom, and without introducing derivatives of the matter 
terms. In Section 4.3 we will introduce an analogous re-formulation of the 
ADM equations. 


4-2 Hyperbolic formulations 


In light of the disadvantages of the ADM system, a large number of 3 + 1 hy¬ 
perbolic formulations of general relativity have been developed recently (see 
the recent review article by |Rcula| (|1998|) for an extensive survey and a more 
complete list of references). The first formulations that departed from the as¬ 
sumption of harmonic coordinates (|Choquct-Bruhati |1952| , |1962|; [Fischer and 
Marsden|. |1972|) were based on a spin-frame formalism ([Fricdrich|, |1981a|,|l5j. 


1985| , |1986a| ,|b|). Other formulations introduced partial derivatives of the met¬ 


ric and other quantities as new independent variables (e.g. [Bona and IVlassd 
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(|1992|) ; |FrittcIIi and Rcula[ (|1994| , |1996|) ; |Bona, Masso, Seidel and Stcla| (|1995| , 
19971 )). In analogy to the electromagnetic example in Section 4.1, Choquet 


Bruhat and Ruggeri| ( |1983|) and [Abrahams et.al .| (|1995| ) took a time derivative 
of equation (2.47) to derive a hyperbolic system that is sometimes referred 
to as the “Einstein-Ricci” system (see also [Abrahams and York| (|1997| )). Al¬ 
ternatively, [Fricdrichl (|1996| ) and |Andcrsom Choquct-Bruhat and York| ( [19971 ) 
used the Bianchi identities to derive another hyperbolic system for general rel¬ 
ativity, sometimes called the “Einstein-Bianchi” system. [Anderson and York 
(|1999| ) developed the “Einstein-Christoffel” system by introducing additional 
“connection” variables. Other recent hyperbolic formulations are based on 
frame or tetrad formalisms (|van Putten and Eardleyl , |1996|; |Estabrook, Robiir 


son and Wahlquist|, |1997|), the Ashtekar formulation (|Yoncda and Shinkai|, 


1999| , p000| , |2001a|; (Shinkai and Yoncda|, |2000|) , a conformal decomposition ( |AL 


cubierre et.al . [ |1999| ), and a so-called A-system that embeds Einstein’s equa¬ 
tions in a larger symmetric hyperbolic system with the constraint surface of 
Einstein’s equations as an attractor of the evolution ([Brodbcck et.al. [ |1999|) . 


It is clearly beyond the scope of this article to review all of these formulations. 
Some of these systems have features that are not very desirable numerically, 
in that they restrict the gauge freedom, introduce extra derivatives of the 
matter variables, or introduce a large number of auxiliary variables. Only few 
of these formulations have been implemented numerically, and most of these 
implementations assumed certain simplifying symmetry conditions (e.g. spher¬ 
ical symmetry). Some of these implementations showed advantages over the 
ADM formalism (e.g. [Bona and Masso| ( |1992|) ), but others also revealed addi¬ 
tional problems. |Schccl et.al. \ (|1997| . |1998| ), for example, found that a particular 
equation in the “Einstein-Ricci” system produced an instability, which could 
be removed in spherical symmetry, but not in more general 3D simulations. 
Very few hyperbolic systems have been implemented in 3D, including that of 
Bona and Masso ([Bona and Masso| , |1992| ; [Bona, Masso, Seidel and Stcla| , |1995| , 
1997]), and generalized versions of the “Einstein-Christoffel” system ([Anderson 


and York] . |1999|) . The latter have been implemented numerically by [Kidder 


Scheel and Teukolskyl (|2001[ see also [Kidder et. al\ (|20QQ|) ) using spectral meth¬ 
ods. Since this formulation is particularly elegant and currently seems like the 
most promising hyperbolic formulation, we provide a brief summary. Some 
properties of this formulation have been analyzed by Calabrese, Lehner and 
Tiglio| (|2002|) ; |Calabrcse et.al . | Q2002|) and [Lindblom and Scheel ( [2002| ) . 


Starting with the ADM formalism as presented in Section 2, and adopting the 
notation of [Kidder et.al .| ([2000D, we define the new variables 


fkij T r [ij] m T 7fcj7 


,lmj^ 


[li]m • 


(4.11) 


These functions are now promoted to independent functions. It can then be 
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shown that the evolution equations (2.46) and (2.47) can be rewritten as 


dt^jij 2aKij 
d t Kij + aj kl difkij = OiMij 
dtfkij 4 ~ (xNfoij . 


(4.12) 


Here we have used the abbreviation 

d,=d,~ Cg 


(4.13) 


and the source terms M tJ and N ijk are given by 

Mij = 7 k \K k iKij - 2 K ki Kij) + 'y kt 'y mn (^fk mi f[in]j 

^^fkm[nfl\ij fikmfjln 4 ” 8 f(ij)kf[ln]m 4 “ 4 fkm(ifj)ln 

8 fklifmnj 4 “ 20 mn ~ 13 fikifj mn) 

-didj In a - (<9. t In d) (dj In d) + 27 ^ 7 ^ 7 ”™ (/ femn <9; In d 
- fkmidn In d) + 7 W ((2/fe)fe - fkij)di lnd 

+4 fki(idj) In d - 3(/ ifcZ <9 ? In d + f m di In a)) 

— 8nSij + ATT'-jijT 


(4.14) 


and 

N kij = 7 mn (4 K k{i fj) mn / ^fmn(iKj)k 4“ ^iji^fmnk 3 f kmn)^ 

+2 7 mn 7 pq ( K K mp ('y k{i f j)qn - 2f qn{ rf j)k ) 

~^’ r yk(i-^j)m{^ l fnpq 6 fpqn) 4 “ pq(ihj)k ^ 7 k(ifj)pq)^ ( 4 . 15 ) 

K-ij^k In d H~ 27 ((iTj) k lnd n 7fc(i ) In d) 

+ 167T7 fc (iJ i ). 


We have also used the “densitized” lapse function 

d = 7 -1 / 2 a (4.16) 

and, in addition to the matter terms defined in Section 2, the four-dimensional 
trace of the stress energy tensor 

T = g ab T ab . (4.17) 
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The above equations form the “Einstein-Christoffel” system of |Anderson and| 
York| (|1999 ). Evolutions of a single black hole using a spectral implementation 


of this system are still unstable (compare Section 8.1), but [Kidder, Schccl and] 
Teukolskyl (|2001[) were able to show that the lifetime of these simulation can be 


extended to late times by a certain generalization of the equations. This gener¬ 
alization involves the redefinition of variables and an addition of constraints, 
and embeds the above equations into a 12 -parameter family of hyperbolic 
formulations. The stability properties of the system depend strongly on the 
choice of the free parameters, which can be understood analytically in terms 
of energy norm arguments (|Lindblom and Scheel|, |2002|). 


The first order symmetric hyperbolic system (4.12) is equivalent to the original 
set of evolution equations (2.46) and (2.47). Since the f ki j are evolved as 
independent functions, their original definition (4.11) can be considered as 
a new constraint equation in addition to (2.44) and (2.45). This system is 
particularly elegant because the source terms M t] and on the right hand 
sides do not contain any derivatives of the fundamental variables (other than 
the arbitrary lapse function a). Equations (4.12) can be combined to yield a 
wave equation for the components of the spatial metric 7 ^ in which the right 
hand sides appear as sources. 


Hyperbolic systems have the great advantage that their characteristic struc¬ 
ture can be analyzed. In the system (4.12), all characteristic fields propagate 
either along the light cone or normal to the spatial foliation. This knowledge 
can be used for the construction of boundary conditions, both at the outer 
boundaries and on inner boundaries if the interior of black holes is excised 
(see, e.g., [Kidder et.al .| (|2000|) and Section 8.1). 


4-3 The BSSN formulation 


Following the electrodynamic example of Section 4.1 we can also eliminate 
the mixed second derivatives in the Ricci tensor with the help of auxiliary 
variables ( |JNakamura, Oohara and Kojima| , |1987|) . In addition, the conformal 
factor and the trace of the extrinsic curvature are evolved separately, which 
follows the philosophy of separating transverse from longitudinal, or radiative 
from non-radiative degrees of freedom (see Section 3). 


We follow here the formulation of [Oaumgarte and Shapiro| (|1999|) , which is 
based on that of Shibata and JNakamura| (1995). We start by writing the con¬ 
formal factor as if = so that 


Tij = e 40 7 ij, 


(4.18) 
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and by choosing it such that the determinant of the conformally related metric 
7 ij is unity, 0 = (ln 7 )/ 12 . As in equation (3.6) we split the trace from the 
extrinsic curvature and conformally rescale the traceless part A,j. Following 
Shibata and Nakamura| (|1995| ) and [Bamngartc and Shapirc| (|1999|) again, we 
choose a conformal rescaling that is different from (3.7) and instead rescale 
Aij like the metric itself 


A i:j = e-^Aij. 


(4.19) 


We will use tildes as opposed to the bars used in Section 3 as a reminder 
of this different rescaling. Indices of A^ will be raised and lowered with the 
conformal metric 7 ^, so that A u = e 4 < M lJ 

Evolution equations for 0 and K can now be found from equation (2.48), 
yielding 

< 9 t 0 = 01 K + /3 l di(j) + \di(3 l (4.20) 

6 6 


and (2.49) 


d t K = -fWjDia + a{AijA^ j + \l< 2 ) + 4tt a(p + S) + ffd t K. (4.21) 

o 

Subtracting these from the evolution equations (2.46) and (2.47) yields the 
traceless evolution equations for 7 ^ 

dtla = -2a Aij + f3 k d k Aij + 7 ikdj(3 k + 7 kj di/3 k - ^^0^. (4.22) 


and A^ 


dtAij = e~ 4 * (-{DiDja) TF + a(B FF - 8nS FF )) 

+a{KAij - 2A il A l j ) (4.23) 

J rf5 k d k Aij + A ik dj(3 k + A k jdi(3 k — ^Aijd k (3 k . 


In the last equation, the superscript TF denotes the trace-free part of a tensor, 
e.g. Rfj F = Rjj fijR/3. Note also that in equations (4.20) through (4.23) the 
shift terms arise from Lie derivatives of the respective variable. The divergence 
of the shift, di/3 l , appears in the Lie derivative because the choice 7 = 1 makes 
0 a tensor density of weight 1/6, and 7 ^ and A tJ tensor densities of weight 
-2/3. 
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According to (3.3) we can split the Ricci tensor into two terms 

Rij = Rjj + Rfj, 


(4.24) 


where Rfj can be found by inserting (j) = In R into (3.3). The conformally 
related Ricci tensor R t] conld be computed by inserting 7 ^ into (2.43), which 
would introduce the mixed second derivatives similar to those in the electro¬ 
dynamics illustration of Section 4.1. In analogy to the new variable T that 
we used to eliminate those mixed derivatives there, we can now define the 
“conformal connection functions” 

r‘ = y*T5» = -7% (4.25) 


where the are the connection coefficients associated with 7 ^, and where 
the last equality holds because 7 = 1 (see Problem 7.7f in |Lightman et.al 
(|1975|) ). In terms of these, the Ricci tensor can be written 


Rij — ~2^ + 7 k(idj)T + Y T(jj)fc + 


^lm 


+ rLr 


k 

innklj 


(4.26) 


The only second derivatives of 7 ^ left over in this operator is the Laplace 
operator 7 * m 7 ijj m ~ all others have been absorbed in first derivatives of P. 
This property of the contraction of the Christoffel symbols has been known 
for a long time ( Pc Dondcr| , |1921| ; Panczos| , |1922|) , and has been used widely to 
write Einstein’s equations in a hyperbolic form (e.g. |LUioquct-Bruhat| (1952); 
Fischer and Marsden| ( |1972| ) ). 


We now promote the P to independent functions, and hence need to derive 
their evolution equation. This can be done in complete analogy to (4.10) by 
permuting a time and space derivative in the definition (4.25) 


a,r‘ = -a, (2 ai« - 27 + -7 


(4.27) 


The divergence of the extrinsic curvature can now be eliminated with the help 
of the momentum constraint (2.45), which yields the evolution equation b j 


d t T l = —2A lj dja + 2a(r i jk A kj — | Y j djK — 8^7^ Sj + 6 A lj dj(f)^ 
AfPdp - Yi<) r r + lY'dj.P + R'AY, + 7 l if3}j. 


(4.28) 


6 Note that the shift terms enter with the wrong sign in equation (24) of Baumgarte 


and Shapiro| (|1999f) . 
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Fig. 4.1. Comparison of the evolution of a small amplitude gravitational wave using 
the ADM equations (dashed line) and the BSSN equations (solid line). The bottom 
panel shows the evolution of the K zz component as a function of time for early times, 
for which both systems agree very well, while the top panel shows the evolution at 
a later time, just before the ADM system crashes. See text for more details. (Figure 
from iBaumgarte and Shapiro| (|1999| ).) 

Equations (4.20) through (4.23) together with (4.28) form a new system of 
evolution equations that is equivalent to (2.46) and (2.47). Since the P are 
evolved as independent functions, their original definition (4.25) serves as a 
new constraint equation, in addition to (2.44) and (2.45). 


While the different formulations are equivalent analytically, the difference in 
performance of numerical implementations is striking. In Figure 4.1 we show 
a particular example from [Baumgarte and Shapiro| (|1999|) . In this example, a 
small amplitude wave ( [Tcukolskyl , |1982|) is evolved with harmonic slicing (see 
Section 5.3), zero shift, and a very simple outgoing wave boundary condition 
(see also |Shiba.ta and Nakamura| ( |1995|) ). Both systems give very similar re¬ 
sults early on, but the ADM system crashes very soon, while the BSSN system 
remains stable. Similar improvements have bee found for many other appli¬ 
cations, including strong field gravitational waves as well as black hole and 
neutron star spacetimes ( |Banmgartc, Hughes and Shapiro| , |1999| ; |Alcubierre 


et. all |2000Io| Jd|; |Lchncr, Huq and Garrisoij , |2000| ; |AIcubicrrc and Briigmannl , 


[2001| ). It is generally found that the ADM system is more accurate initially 
(which is to be expected given that it uses fewer equations), but that the 
BSSN system is much more stable in long term evolution calculations. Given 
this success, the BSSN system, either in the form of |Shibata and Nakamura 
Baumgarte and Shapiro| (|1999|), is currently used very commonly in 


or 


numerical relativity and has been adopted for many recent applications (an in- 

Bata] (|1999a| JE|) ; |Shibata and Uryu| ( |2000a| ) ; |Shibata 


complete list includes 
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Baumgarte and Shapiro ( p() 00 a| .|h|) : |Alcubierre et.al . | ( [2000a| . |2001a| .|b|) ; |Duez 
Banmgartc and Shapiro| ( 2001 ); Font et.al] ( [200 1| )). 


Significant effort has also gone into understanding why implementations of 


the BSSN system are more stable than those of the ADM system flFrittelli 
and Rcula| , |1999| ; |Alcnbierre et.al . | , [2000 b| ; |Millcr| , |2000| ; |Yoncda and Shinkai . 


2002| ; ISarbach et.al] , | 2002 |. see also [Frittclli and Goniezl (|2000|) ). While we are 


still lacking a complete understanding, several arguments point to the propa¬ 
gation of the constraints (compare |Frittclli| ( |1997|) ). [Alcnbierre et.al . | (|2000fo| ) 
linearized the ADM and BSSN equations on a flat Minkowski background and 
showed that modes that violate the momentum constraint propagate with 
speed zero in the ADM equations. They also demonstrate in a model prob¬ 
lem that such zero speed modes lead to instabilities when nonlinear source 
terms are included. Furthermore, they show that by adding the momentum 
constraint in the derivation of the f ® evolution equation (4.28) of the BSSN 
system, the momentum constraint violating modes now propagate with non¬ 
zero speed. Instead of building up locally, as in the ADM system, constraint 
violations can now propagate off the numerical gridQ, presumably stabilizing 
the simulation^. 

This effect can be illustrated very easily with the electrodynamic example 
of Section 4.1 (see [Knapp, Walker and Baumgartc| (|2002| ) ). In the Maxwell 
system, the time derivative of a constraint violation 


C = D l Ei — 47t p e 


(4.29) 


vanishes identically 


d t C = d t (D i E i - 4t rp e ) = D l d t E, - 4vr d t p e 

= -D i D j D j A i + D l D,D 3 A ] - 4tt(DV* + d t p e ) = 0, (4.30) 

where we have used the continuity equation DiJ l +dtp e = 0. Using the modified 
system (4.1), (4.9) and (4.10), on the other hand, it can be shown that C now 
satisfies a wave equation 


d 2 t C = d t D l d t Ei - And 2 p e = <),!)'( -Djl)' Aj + D t T - AnJ t ) - And 2 p e 
= —D i [p j D j d t A i - DAY) - And t (D i J i + d t p e ) 

= D i (D j D 3 (E, + D&) - DiiDWjQ + Anp e )) 


' Assuming that appropriate boundary conditions are imposed. 

8 For simulations of black holes it may be constraint violating modes that prop¬ 
agate along the outward characteristic as opposed to along the normal that cause 
instabilities; see Lindblom and Scheel| ( |2002| ). 
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(4.31) 


= DjD^D'Ei - 4ir p,) = D t D’C. 

The crucial step in achieving this property has been the use of the constraint 
(4.3) to replace D l E i with 47rp e in (4.10). Had we not done that, then the D l E % 
would cancel in (4.31), leading again to a zero propagation speed. Instead, the 
two terms now combine to form C, yielding a wave equation with the speed of 
light as the characteristic speed. 


If C = dtC = 0 initially, then the two systems are equivalent analytically, 
since both will guarantee that C — 0 in the domain of dependence of the ini¬ 
tial surface. Numerically, the two systems behave very differently, since finite 
difference error will lead to a constraint violation \C\ >0. Solving Maxwell’s 
equations in primitive form, such a constraint violation will not propagate 
and will remain constant. As the model problem of |Alcubierre et.al . ] ( |2000b| ) 
suggests, this behavior will lead to instabilities when nonlinear sources are 
included. Using the modified evolution equations (4.1), (4.9) and (4.10), the 
constraint violation C propagates with the speed of light, will leave the numer¬ 
ical grid very quickly, and will ultimately leave behind C = 0. This behavior 
has been verified in numerical implementations of the two systems ( [Knapp 
Walker and Baumgarte| , |2002| ). 


The addition of the constraints to the evolution equations is by no means 
unique to the BSSN system. Pctwcilci] ( |1987| ) pointed out that constraint 
violations can be controlled by adding the constraints to the ADM evolution 
equations. |Frittclli| (|1997|) demonstrated the importance of the propagation 
of constraints in unconstrained evolution calculations, and showed how this 
is linked to adding the Hamiltonian constraint to the evolution equations. 
Other groups have experimented with adding the momentum constraints to 
the ADM equations and have found stabilizing effects ( Petweiler] , |1987| ; |Yoncda| 
and Shinkai , |2001b| ; [Shinkai and Yoncda| , [2002| ; [Kelly etTal . , 200 1|) . Similarly, 
the derivation of many of the hyperbolic systems in Section 4.2 involve an 
addition of the constraints to the original equations. 


5 Choosing Coordinates 


Before the evolution equations as derived in Section 4 can be integrated for a 
set of initial data as constructed in Section 3, suitable coordinate conditions 
have to be chosen. In the framework of the ADM equations, the coordinate 
conditions are imposed with the help of the lapse a and the shift vector f3 l , 
which determine how the coordinates evolve from one spatial slice to the next. 
Often, choices for the lapse are referred to as “slicing conditions”, while choices 
for the shift are called “spatial gauge conditions”. In the following we will 
discuss several different slicings and spatial gauges that have recently been 
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popular in 3 + 1 numerical relativity. Parts of this Section are based on the 
lecture notes of |Baumgartc, Shapiro and Abrahams| (|1998|) . 


5.1 Geodesic Slicing 


Since we are free to choose any lapse and shift, we might be tempted to make 
the particularly simple choice 

a = 1, /T = 0. (5.1) 


For constant lapse the acceleration of normal observers vanishes according 
to equation (2.22), and for zero shift normal observers and coordinate ob¬ 
servers coincide. For this choice of the lapse and shift, coordinate observers are 
therefore freely falling and follow geodesics, which explains the name geodesic 
slicing. 

Example 5.1 In geodesic slicing, the metric (2.60) of Example 2.1 reduces 
to 


ds 2 = — dt 2 + -0 4 (dr 2 + r 2 (d0 2 + sin 2 6d<f> 2 )), 


(5.2) 


which can be identified with the flat Robertson-Walker metric in spherical co¬ 
ordinates. If one further assumes that a homogeneous and isotropic perfect 
fluid is comoving with the coordinate observers, it is easy to show that the 
constraint and evolution equations (2.44) t° (&-4V are equivalent to the well- 
known Friedmann equations. 


Unfortunately, geodesic slicing is a particularly poor choice for most numer¬ 
ical simulations. In a Schwarzschild spacetime every coordinate observer, if 
starting from rest, will fall into the singularity in a finite time, leaving only 
very little time until the numerical simulation would break down (see, e.g., 
the discussion in |5marr and York| (|1978b|) ). To make matters worse, coordi¬ 
nate observers are not only attracted to physical singularities, but also tend to 
form coordinate singularities. This behavior is well known for the evolution of 
even small amplitude gravitational waves. To illustrate this property, imagine 
a small gravitational wave packet located originally at the origin of the coor¬ 
dinate system. After a short time the wave packet will have dispersed, leaving 
behind a flat spacetimes. The coordinate observers will initially be attracted 
by the mass-energy of the gravitational wave packet, and will continue to coast 
toward the origin of the coordinate system even after the gravitational waves 


9 Sometimes these coordinates are also called Gaussian-normal coordinates 
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have disappeared. They will intersect after a finite time and form a coordinate 
singularity. 

This behavior can also be understood from equation (2.49), which reduces to 
d t K = KijK ij = AijA ij + hfi 1 (5.3) 


for geodesic slicing in vacuum. The right hand side is non-negative, and once 
a positive value of K has been induced by any small perturbation it will 
continue to grow without bound. Assuming K = K 0 at t — 0 and A tJ = 0, 
one finds K = 3/i 0 /(3 — K 0 t), indicating that a coordinate singularity forms 
at t = 3/K 0 . Clearly, the usefulness of geodesic slicing is very limited (see also 

; |Anninos et.al .| (|1995cQ ) . 


Shibata and Nakamura (1995) 


5.2 Maximal Slicing 


The above considerations suggest that the pathologies of geodesic slicing can 
be avoided by imposing a condition on the trace of the extrinsic curvature 
K. The most popular such choice is to set K = 0 at all times, which implies 
d t K = 0. Equation (2.49) then yields an equation for the lapse a 

D 2 a = a ( KijK ij + 4vr(p + S)) . (5.4) 

Example 5.2 As we have seen in Section 3.3, the maximal slicing condition 
(5.4) can be combined with the Hamiltonian constraint (3.12) to yield equation 
(3.26) for the product aif. Returning to Example 3.3, we find that all terms 
on the right hand side of (3.26) vanish: A tJ = 0 because of time symmetry and 
equation (3.6), R — 0 because the background metric is flat, and p = S = 0 
because of vacuum. The resulting equation 

D 2 (aif) = 0 (5.5) 


can therefore be solved very easily. Using (3.18) and choosing boundary condi¬ 
tions such that a = 1 for r —» oo and a = 0 on the black hole horizon r = M/2 
we find 


- Mj (2r) 

1 + M/ (2r)' 


With /3 l = 0, the conformal factor given by (3.18) and the lapse by (5.6), the 


39 











metric (2.60) now reduces to 


ds 2 


1 ^ M/M ,,.2 , 

1 + M /(2r) 



(5.7) 


which, as expected, is the well-known Schwarzschild metric in isotropic coor¬ 
dinates. 


It can be shown that maximal slicing extremizes the volume of spatial slices 
spanned by a set of normal observers. A familiar example in Euclidian, three- 
dimensional space is a soap him suspended by a wire loop. Neglecting gravity, 
surface tension will minimize the area of the soap him, and as a consequence 
the trace of its extrinsic curvature vanishes (as Problem 9.31 in [Lightman 
et.al. | (|1975|) demonstrates). 


From equation (2.13) we hnd that in maximal slicing the convergence of normal 
observers vanishes, V a n“ = 0, implying that the normal congruence is conver¬ 
gence free. This property prevents the focusing of normal observers that we 
found in geodesic slicing. 

In strong held regions, this condition will tend to hold back the evolution of the 
slice, and will make proper time “advance more slowly”. This property is called 
singularity avoidance , and is very desirable for many numerical applications 
(see also Section 8.1 and Figure 8.1). Maximal slices starting at v — 0 in a 
Kruskal diagram, for example, never reach the singularity and asymptotically 
approach a limit surface at R — 3M/2 (in Schwarzschild coordinates, see 
[Estabrook et.al. | (|1973|) ; [gmarr and York| (|1978b|) ; [Eardlcy and Smart] ( |1979| )). 
Unlike in geodesic slicing, the entire exterior of the black hole can be covered. 


The properties of maximal slicing have also been studied in a simple model 
problem QSmarr and York| , |1978b| ; |York| . |1979|) . Among other things, this short 
calculation reveals that at late times of an approach to a singularity, maximal 
slicing makes the lapse fall oh exponentially, which is commonly referred to 
as “the collapse of the lapse”. Numerically, this collapse of the lapse has been 
observed in many calculations, including [Evan^ (|1984|) ; |Shapiro and Tcukolsky 
(|1985b| ); [Pctrich, Shapiro and Tcukolsky! (|1985|) . 


While all these properties are very desirable, maximal slicing also has disad¬ 
vantages. Equation (5.4) is an elliptic equation for the lapse function a, and 
in particular in three spatial dimensions solving elliptic equations is very ex¬ 
pensive computationally. Since maximal slicing is only a coordinate condition, 
physical results of a simulation should not be affected if the maximal slicing 
condition is solved only approximately. This suggests that, instead of solving 
an elliptic equation, one could convert that equation into a parabolic equation, 
which is much faster to solve numerically ([Balakrishna et. al. |, |1996|; |Shibata|, 
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1999a|) . As a first step, it may be advantageous to write the maximal slicing 
equation as 


d t K = —cK, 


(5.8) 


instead of dtK , since this prescription will “drive” K back to zero in case 
numerical errors cause it to deviate from zero in the course of some evolution. 
Here c is a positive number of order unity. Inserting (2.49) then yields 

I) 2 a - a (KijK ij + 4t r(p + S)) - f3 i D { K = cK. (5.9) 


This equation may now be solved by rewriting it as a parabolic equation by 
adding a “time” derivative of a 

d x a = D 2 a - a (/\, ; A" + 4vr (p + S )) - (3 i D i K - cK, (5.10) 


where A is an appropriately chosen “time” parameter. If this equation is 
evolved to “late” enough A for every timestep so that convergence d\a = 0 
is achieved, then equation (5.8) is solved and K will be driven toward zero. 
This approximate implementation of maximal slicing is sometimes referred to 
as “K-driver” ( [Balakrishna et.al\ |1996| ) and sometimes as “approximate max¬ 
imal slicing” or AMS (|Shibata| , |1999a ; |Shibata and Uryfl| , |2000a|) . While this 
approach seems very appealing, it also requires some fine-tuning of the pa¬ 
rameters c and A. Moreover, achieving convergence of the parabolic equation 
(5.10) may require a fair number of iteration steps. Most elliptic solvers that 
one would use to solve the elliptic equation (5.4) directly also involve iterative 
algorithms, and in order to save computer time one could limit those to just a 
few iteration steps also. One can therefore think of both approaches as taking 
a few iteration steps toward an approximate maximal slicing solution. 


5.3 Harmonic Coordinates and Variations 

Another particularly simple choice of coordinates is harmonic coordinates , for 
which the coordinates x a are harmonic functions 

VV = 0. (5.11) 

It can be shown that this condition is equivalent to 

( 4) T a = g bc(A) V a bc = 0, (5.12) 
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where is the connection associated with the spacetime metric r/ afe . In¬ 

serting the metric (2.39) into (5.12) shows that in harmonic coordinates the 
lapse and shift satisfy the coupled set of hyperbolic equations 


(d t — ftdj) a = —a 2 K 
(d t - ftdj) ft = - a 2 (fid. In a + y j/ T$ fe ) 


see 


York] Q1979D ). 


(5.13) 

(5.14) 


Harmonic coordinates have played an important role in the mathematical 
development of general relativity, since they bring the four-dimensional Ricci 
tensor d) R ah into a particularly simple form. In the three-dimensional case, 
we found in Section 4.3 that Rij can be written in the form (4.26) in terms of 
the three-dimensional connection functions P = y J '*T l - k . In complete analogy, 
^ R a b can be expressed as 


(4) Rab = -\g cd 9 ab,cd + g C (adb)^r c + ( 4 >r c ( 4 )r (ab)c 

+ 2g*m T c e ^) Vb)cd + 5-^(4) r^( 4 )r ecb 


(5.15) 


Evidently, in harmonic coordinates where = 0 Einstein’s equations reduce 
to a set of nonlinear wave equations ( |L)e Ponder! . |1921| ; |Lanczos| , |1922| ), which 
is why all of the early hyperbolic formulations are based on these coordinates 
(|Choquct-Bruhat| . |1952| ; [Fischer and Marsdcn] . |1972|) . 

Completely harmonic coordinates have been adopted in only few three-di¬ 
mensional simulations ([Landry and Tcukolskyl, |1999| ; |Garfinklc| , |2002|) . More 
popular is so-called harmonic slicing , in which only the time-component ^r 0 
is set to zero. Combining harmonic slicing with zero shift yields a particularly 
simple equation for the lapse, 


d t a = — a 2 Ii , 


(5.16) 


which, after inserting (2.48) for aK , can be integrated to 
a = C^y 1 / 2 . 


(5.17) 


Here C(x l ) is a constant of integration that may depend on the spatial coordi¬ 
nates x l , but not on time. This condition is identical to keeping the densitized 
lapse a = y _1 / 2 a (see equation (4.16)) constant. 

Example 5.3 For the Schwarzschild, metric (5.7) of Example 5.2, Ii, (3 l and 
d t a are all zero. R is therefore obvious that the maximally sliced Schwarzschild 
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spacetime in isotropic coordinates is simultaneously harmonically sliced. For 
numerical purposes, it may be undesirable that the lapse (5.6) vanishes on the 
event horizon at r = M/2. Well-behaved harmonic slices, for which the lapse 
does not vanish on the horizon, have been constructed by \ Bona and Massdj 
\ 1989) and \Cook and Schee\ ( \199Tj ). 


The harmonic slicing condition (5.17) is just about as simple as the geodesic 
slicing condition (5.1), but it provides for a much more stable numerical evolu¬ 
tion (see, e.g., [ghibata and Nakamura| ( |1995|) ; [Baumgartc and Shapiro| ( |1999|) ). 
It does not focus coordinate observers and, as Example 5.3 suggests, met¬ 
ric components of static solutions are independent of time in harmonic slic¬ 
ing, hence allowing for long time evolutions (e.g. |Cook and ScheeT] (|1997| ); 
Baumgartc, Hughes and Shapiro| ( |1999|) ). However, there is no guarantee that 
harmonic slicing will lead to well-behaved coordinates in more general situa¬ 
tions (see, e.g. |Alcubicrrc| (|1997|) ; |Alcubierrc and Masso| (|1998|) ; [Khokhlov and 
l\lovikov| ( ^002|) ) and it has been pointed out that the singularity avoidance 
properties of harmonic slicing are weaker than those of, for example, maximal 
slicing (e.g. [ghibata and l\lakamura| (|1995|); cf. |(birlinkl<j (|2002|)). 


Equation (5.17) is an example of an algebraic coordinate condition, in which 
the lapse can be found algebraically, without having to solve complicated and 
computer intensive differential equations as for maximal slicing. A generaliza¬ 
tion of that condition has been suggested by [Bona, IViasso, Seidel and Stela 


i [T T Fi^ l 


d t a = —a 2 f(a)K, (5.18) 

where /(a) is a positive but otherwise arbitrary function of a. For / = 1 
this condition obviously reduces to harmonic slicing. For / = 0 (and a = 1 
initially), it reduces to geodesic slicing (Section 5.1). Formally, maximal slicing 
(Section 5.2) corresponds to / —> oo. For f = 2/a the condition (5.18) can be 
integrated to yield 

a — 1 + In 7, (5.19) 


where we have used (2.48) and chosen a constant of integration to be unity. 
This quite popular slicing condition is often called “1 + log” slicing (jBernstcm]. 


|1993|; Bona, IViasso, Seidel and Stela, 

1995|), and has been used in various 

applications (including [Anninos et.al. 

'1995c|);|Bona, Masso, Seidel and Walkei) 

(1998); Alcubierre et.al . | (|2001a); 1'ont et.al. (|2001)). As an algebraic slicing 


condition it has the virtue of being extremely simple to implement and fast to 
solve. It has also been found to have stronger singularity avoidance properties 
than harmonic slicing, which can be motivated by the fact that / becomes 
large when a becomes small, so that it probably behaves more like maximal 
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slicing than harmonic slicing. A similar slicing condition has been suggested 
by Shibata and Nakamura| (|1995|) and Oohara, Nakamura and Shibata ( [1997 ) 
with the goal of enhancing the singularity avoidance properties of harmonic 
slicing. 


5-4 Minimal Distortion and Variations 


In Section 3 we found that the conformally related metric 7 \j has five inde¬ 
pendent functions, two of which correspond to true gravitational degrees of 
freedom and three to coordinate freedom. For a stable and accurate numerical 
evolution it is desirable to eliminate purely coordinate-related fluctuations in 
7 ij, which suggests that one may want to construct a gauge condition that min¬ 
imizes the time change of the conformally related metric. This gauge condition 
is called minimal distortion (see Smarr and York ( 1978a Jo), whose derivation 
we will follow closely). 


In Section 3.3 we introduced the time derivative Uij of the conformally related 
metric, 


u 




= 7 1/3 S.(7- 1/3 7«), 


(5.20) 


(equation (3.19)). Since u tJ is traceless, we can decompose it into a transverse- 
traceless and a longitudinal part 

u ij = u Jj r + u ij, (5.21) 


similar to the decomposition of the traceless part of the extrinsic curvature in 
Section 3.2. The divergence of the transverse part vanishes 

D j u^ = 0, (5.22) 


and the longitudinal part can we written as 

< = DiXj + DjXi - lxjD k X k = (LX)* (5.23) 

3 

Since 7 ^ = 7 ~ 1 /, 3 77 - is a vector density of weight —2/3, the right hand side of 
(5.23) can be identified with the Lie derivative of 7 ^ along the vector X 1 , 

Uij = 7 V3 £x7 v- ( 5 - 24 ) 
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Evidently, the longitudinal part can be interpreted as arising from a change of 
coordinates, generated by X 1 . ft represents the coordinate effects in the time 
development, which can therefore be eliminated by choosing uf 3 to vanish. 
This leaves only the transverse part uf 3 T , which implies that the divergence of 
Uij itself must vanish 


D^Uij = 0 . 


(5.25) 


Combining this with equation (3.20) yields 
D j (Lfl)ij = 2 D^aAj) 


(5.26) 


or 

(A L py = 2 A ij Dja + -n~' j I) jl\ + 16 not?, 

o 


(5.27) 


where we have replaced the divergence of A y ' with the momentum constraint. 
This is the minimal distortion condition for the shift vector f3 l . This condition 
can also be derived by minimizing the action 

A = J u ij u k a lk l :,l 'y 1/2 d 3 x (5.28) 


with respect to j3 l . 

It is also useful to express (5.27) in terms of conformally related quantities. 
Using (L/3yi = and DjS y ' = ? ) for any symmetric, 

traceless tensor, we find 

(A L py = ^- 4 ((A L py + {L^Dj ln^ 6 ) (5.29) 

and hence 

(A L py + (Lp) ij Dj In y 6 = 2 A ij Dja + \af ] D 3 K + 16? x^af. (5.30) 

o 


Here we have used the rescaling (4.19) of Section 4.3 for . 


Evidently, the minimal distortion condition (5.30) is fairly involved and quite 
difficult to solve numerically. |Shibata| (|1999c|) therefore suggested simplifying 
the condition. Just as was the case maximal slicing, simplifying is justified 
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because the condition is only a coordinate condition. We may hope that min¬ 
imal distortion has properties that are advantageous for numerical implemen¬ 
tations, but small modifications of the conditions may still lead to similarly 
desirable properties. 


Shibata ( |1999c|) suggested modifying (5.30) in two steps. In the first step, we 
can express A l i in terms of ( L(3) l i as 


A ij 




(5.31) 


(compare equation (3.22)). Assuming Uij = 0 in this expression, and inserting 
into (5.30) we find 

(A L py = 2 aA ij Dj In(m/r 6 ) + -af j DjK + 16tt ip A af. (5.32) 

3 

Alternatively, this expression can be derived by modifying the original condi¬ 
tion (5.25), D l Uij = 0, which is equivalent to D l ('ijj 2 u i] ) = 0, to 

= 0 (5.33) 


(|Shibata| , |1999c|) . Note that the condition (5.32) is identical to the shift condi¬ 
tion (3.28) that we found in the thin sandwich approach with = 0, K = 0 
and conformal flatness (Section 3.3). The assumptions of the thin sandwich 
approach hence leads to the minimal distortion condition for the shift. 


In the second step |Shibata| ( |1999c|) simplifies the operator (A lP) 1 in (5.32) by 
replacing it with the flat vector Laplacian (A| at /3)*, so that in Cartesian coor¬ 
dinates the covariant derivatives can be replaced with partial derivatives. This 
condition is called the “approximate minimal distortion” (AMD) condition, 
and has been used successfully in various applications ( |Shibata| , |1999a| ,|c]; |Slii-| 
|bata and Uryu| , |2000a| ; [Shibata, Banmgartc and Shapiro| , |2Q00a| JB|) . Strategies 
for solving the flat vector Laplacian are discussed in Appendix B. 


In collapse situations, however, |Shibata| (|1999c|) showed that minimal distor¬ 
tion leads to shift vectors that point outwards, leading to coordinate points 
being shifted outward, and hence to a coarser resolution at the center of the 
collapse. This is clearly not desirable, because one would probably want an 
increasingly fine resolution at the center of the collapse. [5hibata| ( |1999c|) ex¬ 
perimented with various ways of adding an artificial radial component to the 
shift, which compensates this effect and improves the numerical performance. 


A gauge condition that is closely related to minimal distortion is based on 
the conformal connection function f* of the BSSN formulation (Section 4.3). 
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These functions conld be set to zero, which would result in “conformal three- 
harmonic” coordinates (compare Section 5.3). Instead, Alcubicrrc and Briig- 
maim| ([ 200 ip suggested setting their time-derivative to zero, 


d t r = 0. (5.34) 

Inserting this into (4.28) leads to the “Gamma freezing” condition 

7 %- + li‘% + Pdf* - fidfi + fr 

= 2 A^djCi — 2a(T l j k A k i — | ^djK — S 3 + 6A lJ dj(/)j. 


The relation to minimal distortion can be seen by inserting (4.25) into the 
condition (5.34), which yields 

djd t f j = dj(il> 4 u ij ) = 0 (5.36) 


(compare with the divergence conditions (5.25) and (5.33)). To simplify the nu¬ 
merical implementation of condition (5.35), [Alcubierre and Brugmann| ( [2001| ) 
converted the elliptic equation into a parabolic one, 

d t ^ = kd t T\ (5.37) 


where k is positive and where (4.28) should be inserted for d t T l . In analogy 
to the “K-driver” of |Balakrishna et.al. \ (|1996|) (see Section 5.2), this condition 
has been called the “Gamma driver” condition. |Alcubierre et.al. \ ( [1001 bp also 
experimented with a hyperbolic version 


d 2 t (3 l =i;- 4 kd t r l -rjd t f3\ 


(5.38) 


where both k and r] are positive constants. 


Another gauge condition that is closely related to minimal distortion is the 
minimal shear or minimal strain condition, in which the time change of the 
physical metric 7 ^ (as opposed to the conformal metric) is minimized in the 
action (5.28). [Brady, Creighton and Thorncj (|1998| ) suggested to combine this 
condition with a similar condition for the lapse to construct “comoving” co¬ 
ordinates that are particularly well suited for simulations of the slow inspiral 
outside of the ISCO (see also |Garfinklc and Gundlach| (|1999|) and |G arlinklc] 

~eEo 7 ] (PODCD). 
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6 Locating Black Hole Horizons 


A black hole is defined as a region of spacetime out of which no null geodesics 
can escape to infinity. The event horizon , or surface of the black hole, is the 
boundary between those events which can emit light rays to infinity and those 
which cannot. More formally, it is defined as the boundary of the causal past of 
future infinity (see, e.g., [Hawking and Ellis] (|1973|) ; |Wald| ([1984])). It is formed 
by those outward-going, future-directed null geodesics which neither escape to 
infinity nor fall toward the center of the black hole. The event horizon is obvi¬ 
ously a gauge-invariant entity, and contains important geometric information 
about the spacetime. Unfortunately, its global properties make it very diffi¬ 
cult to locate, since, in principle, knowledge of the entire future spacetime is 
required to decide whether or not any particular null geodesic will ultimately 
escape to infinity or not. An event horizon can therefore at best be found 
“after the fact”, meaning after an evolution calculation has evolved to some 
stationary state. 


Locating event horizons after the completion of a calculation may be sufficient 
for diagnostic purposes, for analyzing the geometrical and astrophysical results 
of a black hole simulation, but locating black holes in numerical simulations 
is also important for a more technical reason. The spacetime singularities at 
the center of black holes have to be excluded from the numerical grid, since 
they would otherwise spoil the numerical calculation (see also Section 8.1). 
As we have seen in Section 4, spacetime slicings can be chosen so that they 
avoid singularities (for example maximal slicing, Section 4.2). Typically, how¬ 
ever, these slices quickly develop grid pathologies which also cause numerical 
codes to crash. Realizing that the interior of a black hole can never influence 
the exterior suggests an alternative solution, namely “excising” a spacetime 
region just inside the event horizon from the numerical domain (|Unruh| , |1984| ). 
This approach requires at least approximate knowledge of the location of the 
horizon at all times during the evolution, and the construction of the event 
horizon after the fact is therefore not sufficient. 


In practice one therefore locates apparent horizons during the evolution. The 
apparent horizon is defined as the outermost smooth 2-surface, embedded in 
the spatial slices E, whose outgoing null geodesics have zero expansion. As we 
will see, the apparent horizon can be located on each slice E, and is therefore 
a local-in-time concept. The singularity theorems of general relativity (see, 
e.g., |Hawking and Kllis| ( |1973|) ; |Wald| ( |1984|) ) tell us that if an apparent hori¬ 
zon exists on a given time slice, it must be inside a black hole event horizon. 
This makes it safe to excise the interior of an apparent horizon from a nu¬ 
merical domainf 10 ]. Note, however, that the absence of an apparent horizon 


10 


This is not necessarily true in other theories of gravity. In Brans-Dicke theory, 
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does not necessarily imply that no black hole is present. It is possible, for 
example, to construct slicings of the Schwarzschild geometry in which no ap¬ 
parent horizon exists (|Wald and lycr| , |1991|) . The latter clearly demonstrates 
the gauge-dependent nature of the apparent horizon. In numerical simulations, 
one simply hopes that the chosen slices are sufficiently non-pathological, and 
that the apparent horizon is reasonably close to the event horizon. This is the 
case, for example, for stationary situations, for which the apparent horizon 
and the event horizon coincide. 


Recently, the concept of isolated horizons has also been introduced (see, e.g., 
[Ashtckar, Beetle and Fairhurst| ( |1999| ); |Ashtckar et.al. | (|2000|) ), and the first 
implementations in numerical relativity have been reported in Dreyer et.al. 
( 120021 ) - 


6.1 Locating Apparent Horizons 


Example 6.1 In spherical symmetry, the concept of an apparent horizon is 
quite transparent, and can be illustrated very easily. The apparent horizon is 
defined as the boundary of the region of trapped surfaces, wherein the cross- 
sectional area spanned by a beam of outgoing light rays immediately evolves to 
a smaller area. 


For the metric (2.60), an outgoing light ray satisfies 


dr 

dt 



( 6 . 1 ) 


and the areal radius is given by fi> 2 r. The apparent horizon condition can there¬ 
fore be found from the condition that the total time derivative offfir along the 
right ray vanish 


d_ 

dt 


(if 2 r) 


d_ 

dt 




2rfi> 2 


V 


P 


-fit¬ 
'll) 2 r 


ib 1 

+ 2 ra— + a 


0 . 


( 6 . 2 ) 


The expression in brackets can be rewritten in terms of the extrinsic curvature 


for example, apparent horizons may exist outside of event horizons, so that only a 
part of the region inside the apparent horizon can be excised (see Scheel, Shapircj 
and Teukolskyl ( |l995bl ) for a numerical example). 
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(2.62), which yields 

Jj' 

- ijj 2 rK° 0 + 2— r + 1 = 0. (6.3) 

V 

This condition can be evaluated very easily. For the Schwarzschild solution 
(3.18), for example, one readily verifies that the apparent horizon coincides 
with the event horizon at r = m/2. 

Consider a closed smooth hypersurface of E, and call it S. By construction, 
S is spatial and two-dimensional. Let s a be its unit outward pointing normal 
in E. Obviously s a then satisfies s a s a = 1 and s a n a = 0. Similarly to how the 
spacetime metric g ab induces the spatial metric 'y ab on E (see Section 2.1), the 
latter now induces a two-dimensional metric 


TTlab 'Yab S a S b Qab T Tl a Tl b S a S b 


(6.4) 


on S. We now consider the outgoing future-pointing null geodesics whose pro¬ 
jection on E is orthogonal to S. Up to an overall factor, the tangents k a to 
these geodesics can be constructed, on S, from 

k a = s a + n a , (6.5) 


which automatically satisfies k a k a = 0 and m ab k a = 0. We parallel propagate 
k a away from S with the geodesic equation k a X7 a k b = 0. A marginally trapped 
surface (or sometimes called marginally outer-trapped surface) is now defined 
as a surface on which the expansion 0 of the outgoing null geodesics orthogonal 
to S vanishes everywhere^] 

0 = W a k a = 0 . ( 6 . 6 ) 


The outermost such surface is called the apparent horizon. 

For the purposes of numerical relativity, it is useful to rewrite this equation in 
terms of three-dimensional objects. To do so, we insert (6.5) into (6.6), which 
yields 

0 = V a (n a + s a ) = -K + V a s b (6.7) 


where we have used equation (2.13). The divergence of s a can now be rewritten 


11 Note that this condition is equivalent to m ab V a k b = 0. 
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( 6 . 8 ) 


V a s h = g a b V a s b = ( 7 a fc - n a n b )V a s b = As* - n a n b W a s b 
= As* + n a s b W a n b = As* + s i s j K ij + k a s b W a n b . 

The last term vanishes identically, as can be seen by inserting (6.5) again, 

k a s b V a n b = s b k a W a k b - k a s b V a s b = 0. (6.9) 


We can now combine results to find the apparent horizon equation 

0 = As* — K + s l s j Kij = 0, (6.10) 

which can also be written as 

0 = ml 3 (Asj — Kij) = 0. (6-11) 


(cf. |York| (|1989| ) and |Gundlach| (|1998| )). This condition only depends on spatial 
quantities defined within each slice E, which makes it obvious that it can 
constructed “locally in time” on each slice. 


Example 6.2 We can now verify that condition (6.11) leads to (6.3) in spher¬ 
ical symmetry. In this case, s a only has a radial component, s r = A 2 ? an d 
the only non-vanishing components of nA are m 0e = r y ee and rrW’ = yA The 
term rri lJ K tJ is therefore 

m 11 K tj = K e e + K\ = 2K e d , (6.12) 


and rrW D,Sj reduces to 


m^DiSj = ml 3 (sj t i - T^s fc ) 


2 



r + 



(6.13) 


where we have used the connection coefficients (2.61). Inserting the last two 
equations into the condition (6.11) immediately yields (6.3), as expected. 

Various different methods have been employed to locate apparent horizons. 
Most of these methods have in common that they characterize the horizon as 
a level surface of a scalar function, e.g. 

t(x 1 ) = 0. (6.14) 


The unit normal s a can then be expressed as 
S* = \f j DjT 


(6.15) 
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where A is the normalization factor 


A = (7 v DiTDjT) 1/2 . 


(6.16) 


Inserting these into equation (6.11) yields 
0 = rn i] (\D 1 D,jT - Kij). 


(6.17) 


The expansion 0 is therefore a second order differential operator on r. The 
principal part of this operator is the Laplace operator with respect to the 
two-dimensional metric on S (see also the discussion in |GundlacE| (|1998|) ). 


Most authors have chosen r to be of the for 11 if 77 ] 
t(x‘) = rc(x') - h(0,4>), 


(6.18) 


where rc is the coordinate separation between the point x l and some location 
C l inside the r = 0 surface, and where 9 and 0 are polar coordinates centered 
on C\ In fact, most authors assume that C l is the origin of the coordinate 
system, C l = 0. The function h measures the coordinate distance from C l to 
the r = 0 surface in the direction ( 9 ,0). 

Inserting (6.18) into (6.17) yields a second order differential operator on h(9, 0) 
(which, when considered as a differential operator in two dimensions, is ellip¬ 
tic). In spherical symmetry, where h is a constant, the problem reduces to 
solving an algebraic equation, as we have seen in Examples 6.1 and 6.2. In 
axisymmetry h only depends on 9, and equation (6.17) becomes an ordinary 
differential equation which has to be solved with periodic boundary conditions. 
Solutions have been constructed with shooting methods (for example |Cadez 

), spectral methods ( |Epp- 
fork| , |199U| ). Both spectral 
and finite difference methods have also been employed in three-dimensional 
problems, without any simplifying symmetry assumptions. 

Nakamura, Kojima and Oohara (|1984|) (also [Nakamura, Kojirna and Oohara 
(|1985|) ) adopt a spectral method. They expand h in spherical harmonics 

^max l 

h(9,(j)) = J2 °lmYlm(9 , 0) (6.19) 

0 m=—l 


(|1974p ; PBishop| ( |1984|) ; |Shapiro and Tcukolskyl (|1992 
ley], |1977|) and finite difference methods (|L-ook aiicT 


and construct an iterative algorithm to determine the expansion coefficients 


12 This choice restricts the topology of S to S 2 , and S must furthermore be star¬ 
shaped around C l (compare the discussion in Gundlach ( 1998|) ). 
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ai m , which is most easily explained in the notation of |Gnndlach| (|1998|) . Recall 
that the principal part of the nonlinear operator 0 acting on h is a two- 
dimensional Laplacian with respect to . The key idea is introduce a linear 
elliptic operator, which is easier to invert, and to subtract it from the non¬ 
linear one. [Nakamura, Kojima and Oohara| (|1984|) use the flat Laplacian on a 
2-sphere, 


L 2 h = hfio + cot 9h : e + sin 2 9h 




( 6 . 20 ) 


and rewrite the equation 0 = 0 as 
L 2 h — p0 + L 2 h. 


( 6 . 21 ) 


Here the scalar function p is chosen so that the partial derivative hpg cancels on 
the right hand side (see |Gundlach| (|1998|) for a generalization). Since L 2 Yi m = 
—1(1 + 1 )Yi m , we can now multiply both sides with Yj* m and integrate over S 
to find 


— 1(1 + 1) ai m — I Yim(P ® + L 2 h) dVt, 


( 6 . 22 ) 


where dVt = sin 9d6d(j), and where we have used the orthogonality of the spher¬ 
ical harmonics. A priori, this equation is not very helpful, since the right hand 
side has to be evaluated at S, which depends on the very cq m that we would 
like to determine. However, equation (6.22) can now be used to define an it¬ 
eration procedure, in which the integral on the right hand side is evaluated 
using a previous set of guesses af m to determine a new set a] 1 ^ 1 ■ Obviously, 
this algorithm works only for l > 1, and a oo has to be determined indepen¬ 
dently from the integral on the right hand side alone. A similar scheme has 
been implemented by [Kcmball and Bishop| ( |1991| ) . 


A variant of the spectral method, in which the problem of finding a root 
of 0 is reduced to a multi-dimensional minimization problem, was proposed 
by [Libson et.al . | ( |1996a|) and implemented independently by [Baumgarte et.al. \ 
(|1996|) and [Anninos et.al . | (|1998| ). In this approach, @ 2 is integrated over S 


S = / Q 2 da, 


(6.23) 


where da is the proper area element on S. The function h is again expanded as 
in (6. 19)[~ r:r l, so that S becomes a function of the expansion coefficients a im . A 
standard minimization method is then used to vary the cp m until a minimum 

13 Except that these authors expand h in terms of symmetric traceless tensors in- 
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of S has been found. An apparent horizon has been located if S, according to 
some appropriate criterion, is sufficiently close to zero. 


Shibata (|1997|) (also [Shibata and Uryh| (|2000b|) ) uses the same ansatz (6.21) as 
[Nakamura, Kojima and Uohara| (|1984[ ), but employs a hnite difference method 
instead of a spectral method to solve it. In particular, S is covered with a 
hnite difference grid (6p (pj ), on which h(9, <p) is represented as hij. The oper¬ 
ator (6.20) is then hnite differenced in a straight-forward, second-order fashion 


(L 2 h) iJ = 


hi+i,j 2/pj T j hi-ij 


(a ey 


+ cot 9i 


2A 9 


+ 


-2 n hi,j +1 2 hij + hij -1 


sin _i 9. 


(a <py 


(6.24) 


Equation (6.21) can again be solved using an iterative algorithm. On the right 
hand side, the operator (6.24) can be evaluated for a previous set of values 
hfj. On the left hand side, the same operator acts on the new values 
Evaluating (6.21) at all gridpoints (6p <f>j) then yields a coupled set of linear 
equations for the h™^ 1 , which can be solved with standard techniques of matrix 
inversion (e.g. press et all ( |1992| )). 


Thorneburg (|1996|) and |Huq, Choptuik and Matzncr| ( [2000|) have also imple¬ 
mented hnite difference methods to locate apparent horizons. They, however, 
do not use the ansatz (6.21), and instead hnite difference the nonlinear equa¬ 
tion 0 = 0 directly. The resulting nonlinear system of equations is then solved 
with Newton’s method. A similar method has been used by pchncttcr| (|2002|) , 
except that here the tensor helds on the horizon are represented in Cartesian 
coordinate components. 


Yet another method, a curvature how method, was proposed by |Tod| (|1991|) . 
This method is related to solving an elliptic equation by converting it into a 
parabolic equation. During the evolution in an unphysical “time” parameter, 
the solution of the parabolic problem settles down into an equilibrium solution, 
which is the solution to the original elliptic problem. Similarly, [Tod| ( |1991| ) 
proposes to deform a trial surface S according to 


dx l 

~d\ 


= —s“0, 


(6.25) 


where A is the unphysical time parameter. For time-symmetric data with 
Kij = 0, 0 reduces to the trace of the extrinsic curvature of S in E, -DjS*. The 
apparent horizon then satishes DjS* = 0 and is therefore a minimal surface 

stead of spherical harmonics, which is completely equivalent, but more convenient 
in Cartesian coordinates. 
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(compare Problem 9.31 in [Lightman et.al . | (|1975|) ), for which Tod’s method is 
known to converge. For general data, the flow (6.25) is no longer guaranteed to 
converge, but numerical experience shows that in general it still does. Imple¬ 
mentations of Tod’s method have been described in a number of unpublished 
reports (see, e.g. |GundlacIi| ( |1998|) and {Shoemaker, Huq and Matzner| ( |2000| ) 
for references). 


Gundiach ( |1998|) generalizes the flow prescription (6.25) and describes a new 
family of spectral algorithms which include the methods of |Tod| ( |1991| ) and 
[Nakamura, Kojima and Oohara| ( |1984[ ) as special cases for particular values of 
certain parameters. It is suggested that other algorithms in this family may 
combine the robustness of the former with the speed of the latter. 


Alcubierre et.al. (|2000c|) have compared the algorithm of Pundlach| (|1998|) with 
the minimization method of |Anninos et.al . | (|1998|) for a number of different 
test problems, without finding clear advantages of one algorithm over the 
other. The spectral method of |GundlacIi| (|1998| ) is reported to be generally 
much faster than the minimization routine, but the efficiency of the latter 
could probably be improved dramatically. Both |Anninos et.al\ (|1998|) and 


Baumgarte et.al] (|1996| ) employ Powell’s method for the multi-dimensional 


minimization, since it does not use derivatives of the function and is hence 


easy to implement. |Pfciffer, Teukolsky and Gook (1 

1000) report on replacing 

Powell’s method in the code of Baumgarte et.al. 

(19961) with a Davidson- 

Fletcher-Powcll algorithm (see 

Press et.al. 

(1992 

)) and find a significant speed- 


up. 

Shoemaker, Huq and Matzner (|2000| ) implement a finite difference version of 
a flow description with is also a slightly generalized version of (6.25). They 
also introduce a so-called level flow, in which the flow proceeds to a surface of 
constant, but not necessarily zero expansion 0. The apparent horizon is then 
located by constructing a sequence of such surfaces, including the case 0 = 0. 
This approach may have certain advantages for handling situations in which 
multiple apparent horizons are present. 


6.2 Locating Event Horizons 


Event horizons, formally defined as the boundary of the causal past of future 
infinity, are traced out by outgoing light rays that never reach future null 
infinity and never hit the singularity. In principle, therefore, knowledge of the 
entire future evolution of a spacetime is necessary to locate event horizons. 
In practice, however, event horizons can be located fairly accurately after a 
finite evolution time, once a spacetime has settled down to an approximately 
stationary state. 


55 














































































Hughes et.al. (|1994|) constructed an event horizon finder by evolving null 
geodesics 


z'/ ry-tb ,] ^y>C 

a JL ^ (4)-pa 

Hy? + bc lxlx 


o, 


(6.26) 


where A is an affine parameter. In 3 + 1 form, this equation can be rewritten 
as 


= -aat'iip°) 2 + P\iPkP° ~ (6.27) 

— = 7 ^-/?y, (6.28) 

where we have used p l = dx l /d\ and p° = (Y^PiPj) 1 ^ 2 / a (which enforces 
g ab PaPb = 0). 


From each point in a numerically generated spacetime, for which the lapse 
a, the shift /3 l and the spatial metric 7 ^ are known on a numerical grid, 
light rays can then be sent out in many different directions p l . To determine 
whether or not this point is inside an event horizon, [Hughes et.al. (|1994|) use 
the additional knowledge of apparent horizons. Conceptually, if all light rays 
sent out from a point end up inside an apparent horizon (which is always 
located inside an event horizon) the point is inside the event horizon as well. 
If, on the other hand, at least one light ray sent out from the point escapes 
to large separations, the point is not inside an event horizon. I 11 this way, the 
ejection and propagation of light rays from various points in spacetime can 
determine the location of the event horizon. 


An alternative and quite attractive approach was suggested by |Anninos et.al 
( |1995a|) and [Libson et.al \ ( |1996|) . Realizing that the future directed light rays 


diverge away from the event horizon, either toward the interior of the black 
hole or toward future null infinity, they suggest integrating null geodesics 
backwards in time, which converge to the event horizon. In practice, [I Iuglud 
et.al . | ( |1994|) also employ backwards integration of light rays for the same 
reason. This method is particularly efficient if one can identify a “horizon- 
containing” region in which the event horizon is expected to reside. It is then 
sufficient to integrate light rays from this limited region, and they will quickly 
be attracted by the event horizon. 


Anninos et.al. (|1995a|) and [Libson et.al . | (|1996| ) also pointed out that instead 
of integrating individual null geodesics, one may integrate entire null surfaces 
backward in time. Defining such a null surface as a level surface f(t,x l ) = 0, 
one can rewrite the null condition 


g ab d a d b f = 0 


(6.29) 
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to find an evolution equation for / 


dtf 


-g u dif + yj ( g u dif ) 2 - g tt g ij difdjf 


(6.30) 


An event horizon can then be located by evolving two such surfaces defining 
the inner and outer boundary of the horizon-containing region backward in 
time. The two surfaces will converge very quickly and will bracket the event 
horizon. 


Examples of event horizons in numerically generated spacetimes, for example 
for the head-on collision of two black holes, can be found in [Hughes et.al.\ 
( 1994f) ; [Anninos et.al. ( 1995a ); Matzner et.al. ( 1995 ) and Libson et.al. ( 1996 ). 


7 Binary Black Hole Initial Data 


In this Section we will discuss various approaches to solving the initial value 
equations of Section 3 for spacetimes that describe binary black holes in ap¬ 
proximately circular orbit. Such initial data may be used as initial data for 
dynamical simulations of the plunge and merger, as “snapshots” of the evo¬ 
lutionary sequence up to the ISCO, to locate the ISCO, and finally as back¬ 
ground models in quasi-adiabatic evolutionary calculations to simulate the 
late inspiral phase. For a more complete review of initial data for numerical 
relativity see |Cook| (| 2000 |) . 


Under sufficiently restrictive assumptions, constructing initial data contain¬ 
ing multiple black hole is almost trivial. Recall from Example 3.3 that for 
conformally flat ( 7 ^ = rjij), time-symmetric ( K^j = 0 ) vacuum spacetimes 
(p = j l = 0), the momentum constraint (3.13) is solved trivially and the 
Hamiltonian constraint (3.12) reduces to the simple Laplace equation (3.17), 
which, with suitable boundary conditions, has solutions of the form (3.18). 
Given that the Laplace equation (3.17) is linear, we can obviously add several 
solutions (3.18) and fold 




i + E 

a 


m a 
2 r c 


(7.1) 


for an arbitrary number of black holes. Here rc Q = \\x l — C l a \\, and C l a is the 
coordinate location of the a -th black hole. 

However, even for two black holes this construction is not unique. We could 
define a coordinate mapping equivalent to that of Example 3.3 for each one of 
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the black holes, but would find that for each one the presence of the companion 
would destroy the isometry that we found for a single black hole ( (Linguist 


1963|). In other words, we can think of the two throats as connecting 


one 


sheet or universe with two black holes with two separate asymptotically flat 
sheets, each one only containing one black hole. This topology is called a 
“three-sheeted” topology. The isometry can be restored by adding additional 
throats inside each one of the already existing throats, which correspond to 
mirror images of the companion black hole ( |Misner| , |1963| ; [Kulkarni, Shcplcy 
and York, |1983). This “conformal-imaging” approach leads to a two-sheeted 


topology, in which the two throats connect to identical, asymptotically flat 
sheets. It is also possible to consider a one-sheeted, but multiply connected 
topology, in which the two ends of a “wormhole” represent the two black holes 
(|Misner| , |1960| ) . 


The non-uniqueness of these solutions stems from the fact that Einstein’s 
equations determine the local geometry of a spacetime, but do not fix its 
topology. In more physical terms, solutions with different topologies differ 
by their initial gravitational wave content, which is also not determined by 
Einstein’s equations. 


From an astrophysical point of view, the above solutions are not very interest¬ 
ing because they assume time-symmetry with Kij = 0. For the construction 
of binary black holes in binary orbit, we will be interested in black holes with 
finite momenta, and hence non-vanishing extrinsic curvature. That means that 
the momentum constraint is no longer solved as an identity, and that the ex¬ 
trinsic curvature introduces a non-linear term into the Hamiltonian constraint. 
In the following we will discuss several different approaches to constructing 
such binary black hole solutions. As we have explained in Section 1 , we will 
be particularly interested in binary black holes in quasi-circular orbit. 


7.1 The Bowen-York Approach 


In the Bowen-York approach ([Bowen and Yor~k|, |1980|; [Bowen|, |1982|; |York|. 


1989| ), initial data are constructed using the conformal transverse-traceless 


decomposition of Section 3.2. Assuming maximal slicing (K = 0) and confor¬ 
mal flatness ( 7 ^ = rjij), the momentum constraint (3.13) reduces to 


(A ™W)* = 0 


(7.2) 


(compare Appendix B). This equation is solved analytically by 


^CP = 


4 r< 


-(7 P‘ + n‘ c n{,Pj) 


(7.3) 
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where n l c = (x l — C l )/r is the normal vector pointing away from the center 
of the black hole at C l , and where P l is an arbitrary vector. Constructing the 
extrinsic curvature from this solution then yields 

^cp = Jr (PV C + PV C + (,« + ) . (7.4) 

c 

Inserting this into (3.15) shows that P l is the linear momentum of the black 
hole. 

Since the momentum constraint (7.2) is linear, a binary black hole solution 
can now be constructed by superposition of single solutions 

& = + ^ 2 p 2 (7-5) 


The total linear momentum of this solution is P l = P\ + Pj, and from (3.16) 
we find that the total angular momentum is given by 

J i = e ijk C{P* + e ijk CiPl (7.6) 


The extrinsic curvature can then be inserted into the Hamiltonian constraint, 
which for conformal flatness and maximal slicing reduces to 

A fla V = (7.7) 

Here A flat is the flat Laplace operator. 


7.1.1 The Conformal-Imaging Approach 


At this point, the topology of the solution to be constructed has to be decided 
on. In the conformal-imaging approach a two-sheeted topology is assumed, so 
that an isometry holds across the throats. For this isometry to hold, additional 
image terms have to be added to the extrinsic curvature (7.5) before it is 
inserted into the Hamiltonian constraint (7.7) flKulkarni, Shcplcy and York| , 
I983p . 


Cook ( |1991| ) and |Cook etfal\ (|1993| ) constructed solutions to the Hamiltonian 
constraint for two black holes with arbitrary momenta using the conformal- 
imaging approach. In these simulations, the isometry conditions on the throats 
were used as boundary conditions, so that the singularities inside the throats 
could be eliminated from the numerical grid. The computational disadvan¬ 
tage of this method is that boundary conditions have to be imposed on fairly 
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complicated surfaces. In finite difference algorithms, this can be accomplished 
either with bispherical or Cadez coordinates (|Cadez| , |1971|, see also Appendix 
C of |Cook| (|1991|) ), designed such that a constant coordinate surface coincides 
with the throat, or else with fairly complicated algorithms in Cartesian coordi¬ 
nates. Both approaches, together with a spectral method, have been compared 
in |Cook et.aJ\ (|1993|) . 


To construct equal-mass binary black holes in quasi-circular orbit, one may 
assume that P l = P{ = —P 2 , and that PiC 1 = 0, where C l is the separation 
vector C l = C[ — C l 2 . For a given mass of the binary, the only remaining 
free parameters are the separation C = ||C*||, and the momentum P = ||P*||. 
It is intuitively clear, however, that for each separation C there is only one 
momentum P that corresponds to a circular orbit - namely the one that 
satisfies the equivalent of Kepler’s third law. 


General relativity does not admit strictly circular orbits, since the emission 
of gravitational radiation will lead to loss of energy and angular momentum, 
and hence to a shrinking of the orbit. During most of the inspiral phase (see 
Section 1), however, the separation decreases very slowly, on a timescale much 
larger than the orbital period. It is therefore very reasonable to approximate 
the orbit as quasi-circular. 

To determine this quasi-circular orbit, |Cook| ( |1994| ) suggested locating turning 
points of the binding energy Eb along a sequence of constant black hole mass 
Mbh and constant angular momentum J = CP. Restricting ourselves to non¬ 
spinning black holes, the mass of each individual black hole might be identified 
with the irreducible masspH 

/ A \ 1/2 

Mbh = Mrr ~ > (7.8) 


where A is the proper area of the black hole’s event horizon ( jChristodoulorj 
1970| ). In numerical simulations, A is approximated as the area of the apparent 
horizon (see Section 6.1). The binding energy can then be defined as 

Eb = M — 2M B h, (7.9) 


where M is the total ADM mass of the system, given by (3.14). Lastly, we 
can define l as the proper separation between the two horizons, measuring the 
shortest path from one surface to the other. 


11 The mass of a single black hole in the presence of neighboring black holes is 
not unambiguously defined in general relativity, in contrast to the total mass as 
measured at infinity. 
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Following |Cook| ( |1994|) , quasi-circular orbits correspond to turning points 


dE b 

dl 


M m ,J 


0 . 


(7.10) 


A minimum corresponds to a stable quasi-circular orbit, while a maximum 
corresponds to an unstable orbit. The transition from stable to unstable orbits 
defines the innermost stable circular orbit (ISCO), which occurs at the saddle 
point 


d 2 E b 

dl 2 


Mbh J 


0 . 


(7.11) 


For a quasi-circular orbit, the binary’s orbital angular velocity hi as measured 
at infinity can then be determined fromp^ 




dE b 

dJ 


M B n,l 


(7.12) 


(compare [Friedman, Uryu and Shibata| ( |2002[) ; [Baker] ( p002| )). A motivation 
and illustration of this approach can be found in [Baumgartf] (|2001| ). 

At this point, a word of warning is in order. As we have discussed above, 
relativistic binaries emit gravitational radiation, causing them to slowly spiral 
toward each other, and they hence do not follow strictly circular orbits. The 
very concept of an innermost stable circular orbit is therefore somewhat ill 
defined. Also, the minimum in the equilibrium energy identifies the onset of a 
secular instability, while the onset of dynamical instability may be more rele¬ 
vant for the binary inspiral (see, e.g., the discussion in [Lai, Rasio and Shapiro 
and [Lombardi, Rasio and Shapiro| ( |1997|) , where it is shown that the 


two instabilities coincide in irrotational binaries). Moreover, it has been sug¬ 
gested that the passage through the ISCO may proceed quite gradually ( Pri| 
and Thorne] p000| ; [Buananno and Da mom] . pOOOj ), so that a precise definition 
of the ISCO may be less meaningful than the above turning method suggests 
(see also Puez et.~aT\ (|2002|) ). Ultimately, dynamical evolution calculations will 
have to simulate the approach to the ISCO and to investigate these issues. 
For the sake of dealing with a well-defined problem, we will here identify the 
ISCO with the saddlepoint of the equilibrium energy (7.11). 


15 Keeping l fixed in this derivative is equivalent to taking the derivative along an 
evolutionary sequence ( |Cook , 1994 ; Pfeiffer, Teukolsky and Cook , 20001) , since the 
difference only appears at second order. I am grateful to H. Pfeiffer for pointing this 
out. 
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Fig. 7.1. The effective potential E^/n as a function of separation l for various values 
of the angular momentum J/pm (thin lines) as obtained in the Bowen-York confor¬ 
mal-imaging approach. Here p is the reduced mass Mbh/ 2 and m is the sum of the 
black hole masses 2Mbh- The thick line connects a sequence of stable quasi-circular 
orbits, which correspond to minima of the binding energy (7.10). This sequence ends 
at the ISCO, identified by the saddlepoint of the binding energy (7.11). (Figure from 
Cool^ (p9^).) 


Combining the conformal-imaging approach with the turning-point method 
|Cook ( 1994 ) constructed the first models of binary black holes in quasi-circular 
orbit (see Figure 7.1, which also illustrates the construction of circular orbits 
and the ISCO with the turning point method). Numerical values for the non- 
dimensional binding energy E h = E^/p, the orbital angular velocity 0 = mil 
and angular momentum momentum J = J/(pm) at the ISCO as obtained by 
Cook| (|1994|) and other authors are listed in Table 7.1. Here p is the reduced 
mass p = Mbh/ 2 and m is the sum of the black hole masses m = 2Mbh- The 


results of |Cook| (|1994| ) were generalized to spinning black holes by |Pfcitfcr, 
Tcukolsky and Cook| Q2000| ). 


7.1.2 The Puncture Approach 


An alternative to the conformal imagining approach is the puncture approach 
suggested by [Brandt and Brugmann ( 1997 ), in which the singularities in the 
Hamiltonian constraint (7.7) are absorbed in an analytic expression. To do 
this, we write the conformal factor -0 as a sum 


0 = u H— 
a 


(7.13) 
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Reference 

E b 

J 


Schwarzschild 

-0.0572 

3.464 

0.068 

Cool (1994) 

-0.09030 

2.976 

0.172 

Baumgarte (2000) 

-0.092 

2.95 

0.18 

Grandclement, Gourgoulhon and Bonazzola (2002) 

-0.068 

3.36 

0.103 

Damour, Jaranowski and Schafer (2000) 

-0.0668 

3.27 

0.0883 

Table 7.1 

Values for the binding energy Ei ,, the angular velocity Q and the angular momentum 


J at the ISCO as obtained in different approaches, 
with 


1 Ad i Ad 9 

- =-(- . 

a 2r Cl 2r C2 


(7.14) 


Iii the limit of infinite separation, the parameters Adi and Ad 2 approach the 
masses of the individual black holes, but in general they are simply constants. 
Since 1/a is a solution to the (homogeneous) Laplace equation, the Hamilto¬ 
nian constraint (7.7) now reduces to 

A flat u = —6(1 + au)~ 7 , (7.15) 

where we have abbreviated 

b = ^a 7 A i:j A ij . (7.16) 


The beauty of this approach is that the poles at the center of the black holes 
have been absorbed into the analytical terms, and that the corrections u are 
regular everywhere ( [Brandt, and h>riigmann| , |1997[ ). Equation (7.15) can there¬ 
fore be solved everywhere, and there is no longer any need to excise the inte¬ 
rior of the black holes from the computational grid. This eliminates the need 
for complicated boundary conditions on the throats, and allows for straight¬ 
forward solutions on very simple computational domains. Since one no longer 
needs the isometry conditions on the throat, one can construct black hole bi¬ 
naries in a three-sheeted topology, so that no additional mirror terms need to 
be added to the extrinsic curvature (7.5), and instead (7.5) can be inserted 
directly into (7.16). The only added complication of this approach is that the 
location of the throat and horizon is not known a priori, and has to be located 
after the fact. This can be done with the apparent horizon finders discussed 
in Section 6.1. 
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Once a solution u has been found, the ADM energy (3.14) can be determined 
from 


M = ~<f>D i ijj d 2 Si = ~ (f D l (-)d' 2 Si- I A ttat u ci 6 x 


2vr 


2tt / 


a 


2tt 


1 f 

: + J^A.2 “I - — / b (1 H- Qf'u) ^d?x. 

27 T J 


(7.17) 


Baumgarte (|2000|) combined the puncture method of [Brandt and Brtigmann| 
(|1997|) with the turning-point method of |Cook| (|1994|) to construct binary black 
holes in quasi-circular orbit. Numerical values for the ISCO are included in 
Table 7.1. |Baker| ( |2002| ) recently found very similar results by keeping a “bare” 
mass fixed instead of the apparent horizon mass. 


7 .2 The Thin-Sandwich Approach 


As we will see in Section 7.3, the numerical results of Fook ( 1994 ) and Baum- 
|garte| (|2000|) disagree with post-Newtonian values for the ISCO by disturbingly 
large factors, which suggests that it would be useful to explore alternatives to 
the Bowen-York approach. 


This has been done recently by |Gourgoulhom Grandclement and Bonazzola 
(|2002| ) and |Grandclement, Gourgoulhon and Bonazzola] (|2002|) , who have con¬ 
structed binary black holes in the thin-sandwich approach (Section 3.3). In the 
thin-sandwich formalism, equations (3.27) - (3.29) are solved for the confor¬ 
mal factor i/j, the lapse a and the shift (3 1 . Unlike in the Bowen-York approach, 
the extrinsic curvature is not constructed from the analytical solution (7.5), 
but is instead computed from i/j, a and f3 l (equation (3.22)). It is important to 
note that |Gourgoulhon, Grandclement and Bonazzola] ( |2002| ), just like |Gook 
( 1994|) and Baumgartc| ( [20001) , assume maximal slicing and conformal flatness. 


Following |Cook| ( |1991|) , |Gourgoulhon. Grandclement and Bonazzola! ( p002|) 
and |Grandclement, Gourgoulhon and Bonazzola] ( |2002 ) adopt the conformal- 
imaging approach (Section 7.1.1), which means that suitable boundary condi¬ 
tions have to be imposed on if), a and (3 l on the throat. This leads to a technical 
problem. The isometry condition requires that both a = 0 and /3 r = 0 on the 
throat. For a = 0 (and i/j non-zero) in equation (3.22), the extrinsic curvature 
can only be regular if the derivative of /3 r also vanishes. The shift therefore has 
to satisfy more conditions than can be imposed. [Grandclement , Gourgoulhon| 
and Bonazzola) Q2002D circumvent this problem by introducing a regulariza¬ 
tion for the shift which then slightly violates equation (3.28), and proceed by 
hoping that this violation is small (see also the discussion in |Gook| (|2002|), who 
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Conformal factor (Z=0) 


Lapse function (Z=0) 


Shift vector (Z=0) 



x [km] 


x [km] 


x [km] 


Fig. 7.2. Isocontours of the lapse function a, the conformal factor ijj and the shift 
vector /?* in the orbital plane 2 = 0 at the ISCO as obtained with the thin-sandwich 
approach. The thick solid lines denote the surfaces of the throats. The kilometer 
scale of the axis corresponds to an ADM mass of 31.8 Mq. (Figure from Grand- 
clement, Gourgoulhon and Bonazzola| (|2002[ ).) 


introduces an alternative set of boundary conditions). 


'Fhe approach of Grandclement, Gourgoulhon and Bonazzola 

(2UU2 

) differs 

in another way from the calculations of 

Cook 

(1994|) and |Baumgarte| (2000), 


namely in the determination of circular orbits. Instead of adopting the turning- 
point method, they compute the Komar mass ( |Komar| , |1959|) 


M k = f 7 ij (D ia - f3 k K ik )d 2 Sj (7.18) 


in addition to the ADM mass M (3.14). In general, the two mass definitions 
lead to different masses, but they agree in stationary spacetimes ( |Bcig| , |1978| ). 


Since quasi-equilibrium solutions approximate stationary spacetimes, [Grand- 


clement, Gourgoulhon and BonazzoIa| (|2002|) suggest that quasi-circular orbits 


can be identified by demanding 


M = M k . 


(7.19) 


The angular velocity 0 enters the equations (3.27) - (3.29) only through 
the outer boundary condition on the shift flGrandclement, Gourgoulhon and| 
BonazzoIa|, |2002|, equation (9)). This value can be varied until M = Mk has 


been achieved. Note that this condition cannot be applied in the Bowen-York 
approach, since it does not provide the lapse and shift that is needed to eval¬ 
uate the Komar mass (7.18). 


Instead of using equation (7.12) to determine fl, as in the turning-point 
method, it can now be used to construct an evolutionary sequence. The ISCO 
can be identified by locating a simultaneous minimum in the ADM mass M 
and the angular momentum J. Numerical values are included in Table 7.1. 
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Fig. 7.3. Results for the angular momentum J, the binding energy E^, and the 
orbital angular velocity Cl at the ISCO for different approaches (see also Table 7.1). 
The square marks the Schwarzschild test-mass result, the circle the third order 
post-Newtonian results of Darnour, Jaranowski and Schafej ( |200C| ), the triangle the 
numerical results of Grandclement, Gourgoulhon and Bonazzola] (|2002D, an d the four 
and six-pointed stars the results of Cook| (|1994| ) and~ Baumga,rtc| (|200C| ). The top 
label gives the corresponding gravitational wave frequencies for a binary of two 25 
Mg black holes. Compare Figures 2 in Baumgartc ( [2001 ) and 20 
Gourgoulhon and Bonazzola (2002|). 


m 


Grandclement 


7.3 Comparison and Discussion 


In Table 7.1 we list numerical results from the numerical calculations of ICook 


(|1994|) , IBaumgartc] ( p000| ) and prandclement, Gourgoulhon and Bonazzola 
(|2002|) . We also include the third order post-Newtonian results of P amour J 
Jaranowski and Schafcr| (|2000|, assuming u; sta tic = 0), as well as the analytical 
values for a test particle orbiting a Schwarzschild black hole, Et = \J 8/9 — 1 ~ 
—0.0572, J = 2v^3 ~ 3.464, and D = 1 /6 3//2 ~ 0.0680. The same numbers are 
also plotted in Figure 7.3. 


It is very noticable that the results of |Grandclement, Gourgoulhon and Bonaz^ 


zola 

(2002 

) agree quite well with 

Darnour, Jaranowski and Schafe 

(20001) (see 

also 

Blanchet 

(2002); Darnour, Gourgoulhon and Grandcement 

(2002 

)), and 


There is a dis¬ 
turbingly large discrepancy between these two groups, however, of as much as 
a factor of two in the orbital frequency. It is most likely that these differences 
are caused by the different choices in the initial value decomposition. In fact, 
Pfeiffer, Cook and Teukolsky| (|2002|) recently demonstrated that “different de¬ 


compositions generate different physical initial-data sets for seemingly similar 
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choices for the freely specifiable pieces”, and that “the choice of the extrinsic 
curvature is critical”. The different decompositions seem to lead to different 
transverse-traceless components of the extrinsic curvature, so that the result¬ 
ing data represent physically different slices (see also P amour, Gourgoulhon 
and Grandcement| (|2002|) ). 


This conclusion leads to the question which initial value decomposition leads 
to more realistic initial data, describing binary black holes in quasi-circular 
orbit. As of now we do not have any reliable means of knowing. One way to 
compare the different initial data would be to evolve them dynamically and 
test which one is closer to being in equilibrium. While dynamical evolution 
codes that are reliable enough to make such comparisons are only now being 
developed (see Section 8), preliminary results suggest that evolutions of the 
initial data of |Baumgartel (|2000| ) do not lead to circular orbits QAlcubierrc 
and Seidel , 2002 ). The fact that the two very different approaches of Grand- 
clement, Gourgoulhon and Bonazzola] ( j2002| ) and |D amour, Jaranowski and 
Schafcij ( P000| ) lead to quite similar results also points to these initial data be¬ 
ing more realistic] 117 ]. This is in accord with our arguments in Section 3.3 that 
the thin-sandwich decomposition provides a more natural means of construct¬ 
ing equilibrium initial data (see also |Lousto and Pried ( |1998| )). Together, these 
findings suggest that the currently most promising approach to numerically 
constructing binary black holes in quasi-circular orbit is the thin-sandwhich 
formulation. 


The good agreement between |Gook| (|1994| ) and [BaumgartT] ( |2000| ) may not 
be surprising, since the two approaches differ, other than in their numerical 
approach, mostly in the underlying manifold structure (the conformal-imaging 
approach of |Cook| (|1994|) adopts a two-sheeted topology, while the puncture 
method used by Paumgartc] Q2000I ) leads to a three-sheeted topology). The 
small difference in their results reflects the fact that the strength of the imaged 
mirror poles in the conformal-imaging approach is much smaller than the 
strength of the poles themselves (|Misnci], 


Several future developments would be very desirable. As discussed above, the 
approach taken by |Grandclement, Gourgoulhon and BonazzoIa| ( |2002| ) seems 
very promising. The inconsistency introduced by the regularization of the 
shift, which may spuriously affect their results, could be eliminated either by 
adopting the boundary conditions of |Cook| ( |2002| ), or by combining the thin- 
sandwich approach with the puncture method and hence completely eliminat¬ 
ing the need for inner boundary conditions. 


It would also be desirable to construct binary black holes using yet other 


16 


Even though both Grandclement, Gourgoulhon and Bonazzola (|2002| ) and 
Damour, Jaranowski and Schafer ( 200C| ) make some adhoc assumptions, so that 
their agreement could possibly be purely coincidental. 


67 

















































































approaches. [Marronetti and Matzner| (|2000| , cf. |Pfeiffer, Cook and Teukolsky 
(|2002|) ) suggested a construction based on black holes in Kerr-Schild coor¬ 
dinates (see Section 8.1), which would avoid the assumption of conformal 
flatness. This approach has the additional benefit that the slices are regular 
across the horizon, which is very desirable for dynamical simulations. Re¬ 
sults for quasi-circular orbits, however, are not yet available. Another possi¬ 
ble approach that avoids the assumption of conformal flatness is to adopt a 
post-Newtonian solution as a background solution, and to solve the Hamilto- 
nian and momentum constraints for corrections to the post-Newtonian metric 
flTichy et.al.[ |2002| ). 


8 Dynamical Simulations of Binary Black Holes 


The dawn of numerical relativity dates back to the calculation in axisynnuetry 
of the head-on collision from rest at large separation of two identical black 
holes by Larry Smarr and his collaborators (see, e.g., |Smar I] (|1979|) ; |Anninos| 
et.al. | (|19931) )• The result was that the emitted radiation is roughly 0.1 % 
Me 2 , significantly less than the upper limit of 29 % allowed by Hawking’s 
area theorem, but in accord with expectations from strong-field perturbation 
theory. 

While significant progress has been made recently in the dynamical evolution 
of binary black hole scenarios in full 3D, even simulations of single black holes 
are still facing stability problems that are only poorly understood. A likely 
candidate for the origin of these problems is the handling of the black hole 
singularities, which we will discuss in Section 8.1. In spite of these difficulties, 
several groups have performed preliminary simulations of binary black hole 
mergers. We will review these efforts in Section 8.2. 


8.1 Singularity Avoidance and Black Hole Excision 


Simulations of black holes are greatly complicated by the presence of singulari¬ 
ties. Encountering such a singularity during a computation would clearly have 
very unfortunate consequences for the numerical simulation. Two approaches 
have therefore been suggested to avoid these singularities: using singularity 
avoiding coordinates, and excising the black hole interior. 

Traditionally, many black hole simulations, in particular in spherical and ax¬ 
ial symmetry, have used singularity avoiding coordinates to model black hole 
spacetimes. This approach is illustrated in Figure 8.1, in which a family of 
singularity avoiding time slices wrap up around the newly formed singularity 




























Fig. 8.1. A black hole spacetime diagram showing various singularity avoiding time 
slices that wrap up around the singularity inside the horizon. Such slicings allow 
short term success in the numerical evolution of black holes, while at the same time 
causing pathological behavior due to “grid stretching” that eventually dooms the 
calculation at late times. (Figure from |Anninos et. al\ (1995b|).) 


inside the horizon. Among the slicing conditions that avoid singularities is 
maximal slicing, which we discussed in Section 5.2. Another singularity avoid¬ 
ing slicing condition that proved quite useful in spherical and axisymmetric 
simulations is polar slicing ( |Eardlcy| , |1982| ; [Bardeen and Pira~n| , |1983| ; |Shapiro| 
and Tcukolskyl , |1986| ; ISchindcr, Bludman and Piran| , |1988| ). The properties of 
maximal and polar slicing have been analyzed for Oppenhcimer-Snyder col¬ 
lapse in Pctrich, Shapiro and Tcukolskyl (|1985| , |1986| ). In three-dimensional 
simulations, the “1+log” slicing (see Section 5.3) has been used fairly com¬ 
monly. ft has been shown to have singularity avoiding properties not unlike 
maximal slicing ([Bona, IVlasso, Seidel and Stela| , |1995| ; |Anninos et.al . | , |1995c| ), 
and as an algebraic condition it is much easier to implement. 


The problem with singularity avoidance is that the time slices become increas¬ 
ingly pathological, as illustrated in Figure 8.1. For late times, they have to 
“stretch” all the way back to an early time in order to avoid the singularity. 
This “grid stretching” will lead to very steep gradients, which will lead to 
large numerical error, and which will ultimately cause the code to crash. One 
could try to resolve these steep gradients with an increasingly large number 
of gridpoints. This attempt leads to a “grid sucking” effect, and is similarly 
undesirable. Instead of solving the problem it only postpones the code crash 
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to a slightly later time, and moreover leads to large amounts of computational 
resources being spent on uninteresting regions of the spacetime. 

If, unlike in Figure 8.1, the initial time slice already contains a singularity, it 
may be possible to absorb the singular part of the solution into an analytical 
expression, similar to the puncture method of [Brandt and Briigmann| ( [19971) 
(see Section 7.1.2). For a single black hole, for example, one can factor out the 
time-independent conformal factor (3.18), and evolve only the conformally re¬ 
lated metric, which is regular everywhere. Functions and their derivatives can 
then be evaluated from the numerical functions together with the analytical 
factors, as long as the singularity of (3.18) does not happen to lie on a grid 
point. 


In the first dynamical simulation of a black hole in three spatial dimensions, 
[Anninos et.al . | (|1995c|) adopted the traditional ADM formulation (Section 2) 
together with geodesic, “1+log” and maximal slicing (Sections 5.1, 5.2 and 
5.3). The simulations in geodesic slicing encounter the central singularity very 
early on. In maximal and “1+log” slicing grid-stretching effects produce very 
large gradients very rapidly, as expected from simulations in spherical or ax¬ 
ial symmetry, and the code crashes after a fairly short time (up to about 
50M in some cases). [Anninos et.al . | ( |1995c| ) experimented with different ways 
of treating the singularities at the center of the black hole, including impos¬ 
ing isometry conditions (see Section 7.1.1) on the black hole throats and the 
analytical factoring described above (see also [Briigmami| (|1996[) ). 


Unruh (|1984[ ) pointed out that there is no need to include the interior of 
the black hole in the numerical simulation, because by the very definition of 
an event horizon (Section 6.2) the exterior is causally disconnected from the 
interior. If no information can propagate from the interior to the exterior, 
there is no need to numerically simulate the interior^} From Figure 8.1 it is 
evident that black hole excision may avoid both the singularity at the center of 
the black hole as well as the grid stretching effects associated with singularity 
avoiding timeslices. 


Black hole excision is currently considered the most promising approach to 
dynamical black hole simulations. It has been successfully implemented in 
spherical symmetry (e.g. [Seidel and Sucn| (|1992|) ; |S dieel. Shapiro and Teukol 


and [Lclincr et.al . | ( |2000|) ), but most implementations in three 
spatial dimensions are still plagued by instabilities. 


Various design choices have to be made in the implementation of black hole 
excision. As discussed in Section 6, the location of the event horizon is un¬ 
known during the dynamical simulation, so the excision is instead based on 


17 Unless, of course, one is interested in the structure of the singularity; see, 
Berger and Moncrief ( 1993| ); Hiibnci ( 1996 ); Garfinkle| ( [2002| ). 


e.g. 
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the location of the apparent horizon. Typically not all gridpoints inside the 
apparent horizon are excised, and instead a buffer zone of a few gridpoints 
is left inside the black hole. At each timestep, the grid points on the inner 
boundary have to be filled with numerical data. If a hyperbolic formulation of 
Einstein’s equation is used with characteristics that lie either on or inside the 
light cone (e.g. the Einstein-Christoffel system discussed in Section 3.2), then 
all these characteristics must be outgoing on this inner boundary (i.e. toward 
the singularity), as is evident from the tilting of the light cones in Figure 8.1. 
In this case, the evolution equations can be used to fill the inner boundary 
points (“boundary without a boundary condition”). 


The tilting of the light cones points to another computational problem. While 
a normal observer n a always lies inside the light cone, a coordinate observer 
t a may lie outside the line cone, resulting in “super-luminal” grid velocities. 
This happens when the shift vector /3 l becomes large (see equation 2.31). As 
a consequence, a gridpoint x Tj) 1 on a timeslice t n+l may be causally discon¬ 
nected from the same gridpoint x^ k on the previous timeslice t n (see Figure 1 
in [Scheel et.al .| (|1997| ) for an illustration). Some numerical schemes become un¬ 
conditionally unstable in such a situation, while others are stable only for very 
small time steps At. To avoid this problem, several groups have used causal 


differencing ([Seidel and Suen, 

1992|; [Alcubierre and Scliutz, 1994]; IScheel et.al, 

|l997|; |Cook et.al . |, [1998 

; Gundlach and Walker, 

1999|; [Lchner, Huq and Garri- 


the gravitational fields along n a instead of f a , and then use an interpolation 
to shift the fields along the shift f3 l . 


If spatial derivatives have to be evaluated in the updating scheme, then deriva¬ 
tives in directions orthogonal to the surface of the excised region have to be 
computed using extrapolation or one-sided schemes on the inner boundary. 
In general, the coordinates will not be aligned with this surface, so that an 
appropriate algorithm has to be constructed. 


If a black hole moves through the computational grid, then grid points re- 
emerge from the excised region on the trailing side of the black hole. These 
grid points have to be filled with data points, presumably by extrapolation. 




) implemented black hole excision together with causal dif- 
i static and boosted black holes, using the ADM formalism, 
hey adopted Schwarzschild black holes in Kerr-Schild form 

Cook et. 
ferencinj 
As initia 
(compar 

al. ( 
for 
1 da 

1998 
botl 
ta, t 

e [Jhandrasekhar (1992); 

Marsa and Choptuikj (|1996|); iVfatzner, Huq 

and Shoemaker ( 

1999); 

Marronetti and Matzner (2000|); |Yo, Baumgarte and 

Shapiro 

(2DDTB)). 


In ingoing Kerr-Schild form, the Kerr metric is given by 

ds 1 = g a bdx a dx b = ( 7] ab + 2 Hl a l b )dx a dx b , (8.1) 


71 






















































































Fig. 8.2. Metric component g zz along the 2 -axis as a function of time. The flat region 
that moves diagonally to the right represents the excised region (inside the black 
hole). Note that points at the trailing edge (left side) are smoothly updated as the 
hole moves toward positive 2 . Coordinate effects are seen to appear near the inner 
boundary. (Figure from Cook et.al. ( 1998|) .) 


where H = H(x a ) is a scalar function. The vector l a is null both with respect 

to flab and CJabi 


>A,i 6 = g^ljb = 0 . 


( 8 . 2 ) 


Comparing with the ADM metric (2.41), we can identify 


a = (l + 2Hl 2 t )- 1/2 
A — 2 HltU 

r Yij Cj + 

For the time-independent Schwarzschild spacetime we have 


(8.3) 

(8.4) 

(8.5) 


H = M/r 

l a =(l,Xi/r), 


( 8 . 6 ) 

(8.7) 


where M is the total mass-energy and r 2 = x\ + x 2 + x 3 . This form is equiva¬ 
lent to the ingoing Eddington-Finkelstein form (compare [Misner 


Thorne and 


W heelerj (|1973|) ). The Kerr-Schild form (8.1) is invariant under boosts, and 
therefore can be used to represent either a static or boosted Schwarzschild 
(or Kerr) black hole. The other great advantage of Kerr-Schild coordinates 
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Fig. 8.3. The root means square of the change of the lapse function a from one time 
step to the next, in the simulations of Alcubierre and Briigmann (2001) for “1+log” 
slicing and static shift. Shown are the results for two runs, one performed in octant 
symmetry, and an identical run except that it does not assume any symmetry. The 
former evolves stably until very late times, while the latter develops an instability 
after a few 100M. (Figure from Alcubierre and Briigmann (2001).) 


is that they are regular on the horizon and hence extend smoothly into the 
black hole. In Schwarzschild coordinates, for example, the lapse vanishes on 
the horizon. As a consequence, infalling radiation will never cross the horizon 
in Schwarzschild coordinates, and will instead “pile up” just outside the hori¬ 
zon. This has unfortunate consequences for numerical simulations, because 
this piling up will produce increasingly small length-scales which cannot be 
resolved with any given finite-difference resolution (see [Rezzolla et.al . | ( j!998|) 
for an illustration in three spatial dimensions). 


Rezzolla et.al. (1998 


Cook et.al. (|1998| ) use analytic Kerr-Schild data as initial data, and adopt 
the Kerr-Schild lapse and shift to fix the coordinates. For static black holes, 
evolutions of up to 95M are achieved, a slight improvement over the results of 
Anninos et.al . | ( |1995c|) (compare also with the run-time of 138 M achieved by 
33)- More importantly, |Cook et.al. ] (|1998| ) demonstrate that black 
advected through a numerical grid. Figure 8.2 shows the excised 
region inside the black hole moving through the grid, and also demonstrates 
that points at the trailing edge, which re-emerge from the black hole, are 
updated smoothly. 


Daucs| (|1996|) | 
holes can be 


18 Note also that characteristic evolutions of black holes in three spatial dimen¬ 
sion had already achieved runtimes of about 1400M flGomez et.al.[ |199S| ). However, 
modeling binary black holes with characteristic methods requires multiple match¬ 
ing between characteristic and presumably Cauchy methods, which has not been 
established yet. 
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A new and quite successful approach was suggested by |AIcubicrrc and Briig-| 
mann| Q2001 , see also Alcubierre et.al. ( |2001b )). They use the BSSN formalism 


(Section 3.3) instead of the ADM formalism, and a particular simple black hole 
excision scheme. Instead of excising a (nearly) spherical region, they excise a 
cube, which is well adapted to the Cartesian coordinates. Instead of trying to 
construct boundary conditions from the evolution equations, they use a simple 
but stable boundary condition on the inner boundary. Finally, instead of using 
causal differencing, they use centered differences except for advection terms 
on the shift (i.e. terms involving /T<9j), which are differenced with an upwind 
scheme along the shift direction. A modification of this scheme was suggested 
by |Yo, Baumgarte and Shapiro] ( |2001b|) . 


Alcubierre and Briigmann (|2001|) use a static black hole in Kerr-Schild coor¬ 
dinates as initial data, and experiment with several slicing and gauge condi¬ 
tions, including using the analytical lapse and shift of the Kerr-Schild solution, 
maximal slicing, “1+log” slicing, and “Gamma-freezing” (see Section 5.4). An 
example of their results is shown in Figure 8.3, where the root mean square 
of the change of the lapse from one time step to the next is plotted as a 
function of time for “1+log” slicing and static (analytic) shift. Two findings 
are quite remarkable. When evolved in octant symmetry the evolution settles 
down exponentially until the changes in the lapse reach machine precision. 
No instability is encountered. However, when octant symmetry is relaxed, an 
otherwise identical evolution develops an instability after a few hundred d/” . 
Similar results were reported by |Schccl| ( [2000| ) , who use a completely indepen¬ 
dent implementation (namely spectral methods) of a different formulation of 
Einstein’s equations (namely the Einstein-Christoffel system described in Sec¬ 


tion 3.2). Improvements over these results have been reported by [Laguna and 
Shoemaker] (|2002| ); [Lindblom and Scheel| (|2002| ); |Yo, Baumgarte and Shapirc 


(|2002|) ; |Alcubierre et.al . | (|2002|) . In the BSSN calculation of [To, Baumgarte 
|and Shapiro| ( |2002| ), the octant symmetry has been removed and long-term 
stability has been achieved both for static and rotating black holes. 


8.2 Evolution of Binary Black Holes 


Disregarding the stability problems in simulations of single black holes, some 
groups QBriigmanri |1999| ; |Brandt et.al . | , [2000| ; |Alcubierre et.al . |, |2001a| ; |Baker 
et.al\, [ 200 1 |, | 2002 |) have initiated preliminary simulations of binary black holes. 


Briigmann (|1999|) , |Brandt eh~a! 1 ( |2000|) and [Alcubierre et.al . | ( |2001a|) simulated 
“grazing” collisions of black holes. In these simulations, two black holes, start- 


19 It is also interesting that when evolved with the ADM formalism instead of the 
BSSN formalism, the simulation terminates after about 30 M , even in octant sym¬ 
metry 
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Fig. 8.4. The merger of apparent horizons in the simulations of Alcubierre et.al. 


( 2001a| ). Shown are marginally trapped surfaces at times 2.5 M, 3.7 M, 5.0 M, and 
6.2 M. The apparent horizon is the outermost of these surfaces. (Figure from Alcu¬ 
bierre et.al\ (2001 a| ).) 


ing out well within the ISCO, are boosted toward each other with a certain 
impact parameter. These initial data do not correspond to the ISCO or any 
other astrophysically motivated configuration, and are instead chosen so that 
a common apparent horizon forms soon after the start of the evolution. 


Brugmann ( |1999|) uses the puncture method of [Brandt and Briigmann| ( [19971) 
both for the construction of initial data (as described in Section 7.1.2) and to 
remove the singularities during the evolution. In analogy to factoring out the 
time-independent term (3.18) for single black holes, |Brugmann| ( |1999|) factors 
out a time-independent term similar to (7.13) for the binary. The remaining 
terms are regular, and can be integrated everywhere without black hole exci¬ 
sion. In this approach the singularities or “punctures”, corresponding to the 
asymptotically flat region inside the black holes, remain at constant coordi¬ 
nate locations during the evolution. [Brugmann] (|1999|) evolves the remaining 
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regular terms using the ADM formalism with maximal slicing and zero shift. 
The simulation terminates very early (at about t = 7 M), and the results are 
too crude to be of any astrophysical interest. They nevertheless show some in¬ 
teresting features, including the formation of a joint apparent horizon, and are 
also noteworthy as the first attempt to simulate binary black hole coalescence 
in three spatial dimensions. 


Brandt et.al. ( f2000|) adopt as initial data a superposition of Kerr-Schild so¬ 
lutions, which serves as an approximate solution to the constraint equations 
(|Marronetti et.al . | , P000| ) . This approximate solution also provides a lapse and 
shift which is used to fix the coordinates during the evolution. The gravita¬ 
tional fields are evolved with the ADM formalism, and the interior of the black 
holes is excised. The results seem to be strongly affected by the outer bound¬ 
aries, which is implemented as a so-called “blending” to the semi-analytic so¬ 
lution. The code again crashes too early (after about t = 15 M on the “large” 
computational domain) to provide any interesting astrophysical results. 


Alcubierre et.al. (|2001a|) return to the approach of [Briigmann| ( |1999|) , but use 
the BSSN formalism instead of the ADM formalism and the “l+log” slic¬ 
ing instead of maximal slicing during the evolution. The initial data and the 
“analytical” handling of the black hole singularities are very similar to |Briig-| 
|mann| (|1999|) . |Al cu bierre et.al] (|2001a|) take advantage of large computational 
resources at the National Center for Supercomputer Applications (NCSA), al¬ 
lowing them to simultaneously use fairly fine resolution and far outer bound¬ 
aries on a uniform Cartesian grid. With this “combination of resolution, outer 
boundary location and treatment, coordinate choice, evolution system and 
puncture method for the black holes”, [Alcubierre et.al . | (|2001a|) were able to 
achieve evolutions well past t = 30 M (see Figure 8.4 for a visualization of the 
merger of apparent horizons). This is a significant improvement over previous 
results, since this number forms a rough minimum for a meaningful grav¬ 
itational wave analysis. The lowest quasi-normal mode for a Schwarzschild 
black hole has a period of about 17 M, which sets the approximate scale for 
the expected gravitational wavelengths in the ringdown phase. A simulation 
past t = 30 M therefore allows one to capture about two periods of this ring- 
down phase, which [Alcubierre et.al . | ( |2001a| ) have been able to identify. At 
later times, however, grid-stretching effects as well as error originating from 
the outer boundaries start to degrade the numerical results. [Alcubierre et.al. 


(|2002|) have recently reported finding singularity avoiding lapse and shift con¬ 
ditions which permit long-term stable integrations in octant symmetry. 


Progress in the evolution of single and binary black holes has been very rapid in 
the past years (further improvements have already been reported by |Alcubierrc| 
|and Scidc| (|2Q02| ) ), and it is quite likely that this held will continue to develop 
quite quickly in the near future. 
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As an alternative to a fully self-consistent numerical simulation of the entire 
plunge and merger of binary black holes, Baker et.al. ( |2001| . 2002 ) recently 
demonstrated that it may be possible to match to a perturbative “close limit” 
treatment fairly early (the “Lazarus” project; see also |Pricc and Pullin| (|1994|) )- 
Baker etml\ ( |2002|) adopt puncture initial data describing binary black hole at 
the ISCO (see IBaumgartc] ( pOOClj ) and Section 7.1.2), and evolve these with the 
ADM formalism, maximal slicing and zero shift (similar to |Briigmann| (|1999|) ). 
The simulation terminates after about t — 15 — 20M. As an important self- 
consistency test, [Baker etmT. K B ) show that the emitted gravitational wave 
signal and energy is independent, within a certain regime, of the time at 
which the matching is performed. It is likely that this or a similar approach 
will ultimately produce the desired late-time plunge and merger gravitational 
waveforms. 


9 Binary Neutron Star Initial Data 


In this Section we will discuss initial data describing binary neutron stars in 
quasi-equilibrium, quasi-circular orbit. Like the binary black hole solutions 
of Section 7 these models may serve as initial data in dynamical evolution 
equations (see Section 10), for the construction of evolutionary sequences up 
to the ISCO, and as background data for quasi-adiabatic approximations, as 
discussed in Section 10.4. We first discuss hydrostatic quasi-equilibrium in 
general (Section 9.1), and then specialize to corotational (Section 9.2) and 
irrotational binaries (Section 9.3). 


9.1 Hydrostatic Quasi-Equilibrium 


As we have discussed in Section 3.3, it is natural to adopt the thin-sandwich 
formalism for the construction of quasi-equilibrium initial data, since it allows 
us to set explicitly the time derivative of the conformally related metric to zero. 
With very few exceptions (|Usui, Uryit and Eriguchi| , |1999| ; [Usui and Erigucln , 
2002| ), most binary neutron star initial data have also been constructed under 
the assumption of spatial conformal flatness. This has the great advantage 
of dramatically simplifying the equations, and has been justified naively by 
arguing that it approximately minimizes the gravitational wave content in the 
spatial slice (see also the discussion in Appendix C). Under these assumptions, 
the metric is written in the conformally flat form 


ds 2 = —a 2 dt 2 + ifj A r]ij(dx l + /Tdt)(cLc J + ftdt), 


(9.1) 
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and Einstein’s equations reduce to (3.27) through (3.29) for the lapse a, the 
conformal factor if} and the shift f3 l E3 

To determine the matter sources in (3.27) to (3.29) we now have to specify the 
energy-momentum tensor T ab , which, adopting a perfect fluid, can be written 
as 


T ab = Oo +Pi + P)u"u b + Pg 


(9.2) 


Here u a is the four-velocity of the fluid, and po , pi and P are the rest-mass den¬ 
sity, internal energy density and pressure as observed by a comoving observer 
u a . It is also useful to introduce the specific enthalpy 


h = exp 


dP \ 
PO + Pi + P J 


(9.3) 


For a polytropic equation of state 


P 


1+1 In 

KPo 


J 


(9.4) 


where k is the polytropic constant and n the polytropic index[_], the enthalpy 
becomes 

h = Po + Pi + P . (9.5) 

Po 


Note that equations (9.3) and (9.5) imply 
dh = Pq 1 dP. 


(9.6) 


The equations of motion 


Vj ,T ab = 


(9.7) 


yield the Euler equation, which can be written as 
v b V b (hu a ) + V a h = 0. 


(9.8) 


20 


This formalism was developed independently by Wilson and Mathews ( |1989| . 


1995) as an approximate approach for dynamical simulations of binary neutron 


stars. We will discuss this approach in Section 10.2. 

21 It is also common to use the polytropic exponent T = 1 + 1/n. 
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Rest-mass conservation yields the continuity equation, 


V a (p 0 «“) = 0. 


(9.9) 


Following |Shibata| (|199S|), we now write the fluid four-velocity u a as 


u a = u t {l a + V a ), 


(9.10) 


where l a is timelike and normalized so that l* — 1, and V a purely spatial, 
n a V a = 0. We will later assume that l a is an approximate Killing vector. The 
spatial components of the velocity vector in a coordinate system comoving 
with l a , which we will refer to as a corotating coordinate system, then reduce 
to aV 1 . 


The four-dimensional equations (9.8) and (9.9) can now be projected into the 
slice £. Expressing the covariant derivative l b VbU a in terms of a Lie derivative 
along l a and using equation (2.24) yields 


7j“£i(/m a ) + A + UjV J j + V 3 (DjUi - DiUj) = 0 


(9.11) 


for the projection of (9.8) and 

a (C^PqU 1 ) + p Q u t V a l a ') + D i (apou t V' 1 ) = 0 (9.12) 


for (9.9) (see |Shibata| ( |1998| ) and compare footnote 3 in |Gourgoulhon| ( |1998| )) 
Here we have introduced Ui as the spatial projection of hui , 


u.i = 'y i a hu a . 


(9.13) 


So far we have not made any assumptions about the symmetry of the space- 
time. Now, however, we assume that l a is a Killing vector, V Q 4 + Vf,/ a = 0, in 
which case the quantities C\(hu a ), C\{pqU 1 ) and VJ a all vanish in the above 
equations. The Euler equation (9.11) then reduces to the Bernoulli equation 


A f 77 + UjV 3 j + V 3 (DjUi - DiUj ) = 0, 


u 


(9.14) 


and the continuity equation (9.9) becomes 
A(«PoM*W) = 0. 


(9.15) 
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Example 9.1 We can now generalize the spherically symmetric, static solu¬ 
tions of Examples 3.3 and 5.2 to include matter. For a static star we have 
u a = n a so that uf = a -1 . From equation (2.26) we find p = Po + Pi and 
from (2.35) Sy = P^ij and hence S = 3 P. The Hamiltonian constraint (3.12) 
becomes 

D 2 i(> = —27 nf 5 p, (9.16) 


which, in spherical symmetry, can be integrated once to yield 


dip 

dr 



m(r ) 
2 r 2 


(9.17) 


where m{r ) = 47r ffip 5 pr ,2 dr'. In the exterior of the matter source, m(r) is 
independent of r and (9.17) can be integrated once more to yield (3.18). The 
maximal slicing condition (3.26) can be similarly integrated to yield 

d(aip) 27 : r -. . (9 , , amir) 

—Qfi~ = 72 J ^ (P + 6P ) r dr = ( 9 ‘ 18 ) 

with m{r) = 47tck _1 ff otif 5 (p + 6P)r ,2 dr'. For static solutions we have V 1 = 0 ; 
so that the continuity equation (9.15) is satisfied identically and the Bernoulli 
equation (9.14) reduces to 

Dfiah) = 0. (9.19) 


This equation can now be rewritten as 


aDih 


—hDia 



h 

ip 


(Dfiaip) - aDiip) 


or, with (9.6), (9.17) and (9.18), 

dP p$h m + ifi 

dr -0 2 r 2 


(9.20) 


(9.21) 


This is the equivalent of the Oppenheimer- Voikoff equation ( \Oppenheimer ancj 
Volkoff, \193ty) in isotropic coordinates. In the Newtonian limit this equation re¬ 


duces to the Newtonian equation of hydrostatic equilibrium. We point out, how¬ 
ever, that solving the equations of hydrostatic equilibrium for spherical stars is 
far more straight-forward in Schwarzschild areal coordinates then in isotropic 
coordinates presented here (e.g. Misner, Thorne and Wheelei\ <\1973)). 
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For the construction of binaries in quasi-circularp 7 ] orbit we assume that the 
binary has a constant orbital velocity fh The time vector l n in a frame coro¬ 
tating with the binary is therefore a Killing vector and can be constructed 
from the time vector in the binary’s rest frame t a as 

r = t a + nc = an a + p a + nc- ( 9 . 22 ) 


Here we have used (2.31) and have introduced £ a as the generator of rotations 
about the rotation axis. In Cartesian coordinates, £ a = (0, —y, x, 0) represents 
a rotation around the z-axis. The Killing vector l a is sometimes called the 
helicoidal Killing vector ( [Bonazzola, Gourgoulhon and MariT| . |1997|) . 


The Bernoulli equation (9.14) and continuity equation (9.15) simplify for the 
two cases V a = 0 and Djiti — DiUj = 0. The former corresponds to corotating 
binaries, and the latter to the more realistic case of irrotational binaries. We 
will discuss both cases separately in Sections 9.2 and 9.3. 


9.2 Corotational Binaries 


The easier case follows from the assumption that the fluid flow vanishes in the 
frame corotating with the binary 

V a = 0, (9.23) 


which leads to models of corotating binaries. With (9.23), the continuity equa¬ 
tion (9.15) is solved identically, and the Bernoulli equation (9.14) reduces to 

h 

— = const = C (9.24) 

tr 


(see, e.g., Problem 16.17 in |Lightman et.al 


019751 )). 


Using the normalization u a u a = —1, the u t component can be expressed in 
terms of the spatial components tq as 

au l — (l+Y^UiUj^j ^ . (9.25) 

22 Like binary black holes (Section 7.1), binary neutron stars in general relativity are 
not in strict equilibrium. Due to the emission of gravitational radiation, they loose 
energy and slowly spiral inward. For binaries sufficiently far from the ISCO, however, 
the inspiral timescale is much longer than the orbital timescale. This justifies the 
notion of binary neutron stars being in quasi-equilibrium, quasi-circular orbits, as 
we are assuming here. 
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Fig. 9.1. Rest density contours in the equatorial plane for a n = 1 neutron star 
binary close to the ISCO. Each star has a rest mass of M = 0.169, corresponding to 
a compaction in isolation of (. M/R ) oc = 0.175. The contours show isosurfaces of the 
rest-density in decreasing factors of 0.556. (Figure from Baumgarte et.al. (|1998b|) .) 


With V a = 0, and assuming rotation about the £-axis, the four-velocity (9.10) 
reduces to 


u a = «*(!, —fly, fir, 0). 


(9.26) 


Since we are also assuming spatial conformal flatness, (9.25) can then be 
rewritten as 


cm = 


b - ^ ((Sly - n 2 + (Six + 3 y ) 2 + (/3*) 2 )1 


- 1/2 


(9.27) 


Inserting this into (9.24) shows that the enthalpy h is given as an algebraic 
function of a, if) and /3 l , which are determined by the thin-sandwich equations 
(3.27) to (3.29), and the constants f2 and C. From h the fluid variables p 0 , pi 
and P can be computed. Finally, the matter-sources p, j l and S in (3.27) to 
(3.29) can be computed in terms of the fluid variables and the velocity from 
their definitions (2.26), (2.29) and (2.35). 

Self-consistent solutions to the Bernoulli equation (9.24) together with the field 
equations (3.27) to (3.29) can be constructed with iterative algorithms. The 
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Pc 


Fig. 9.2. Rest mass Mq versus central density p c for separations z A = 0.3 (bot¬ 
tom solid line), 0.2, 0.1 and 0.0 (top line) for corotational, n = 1, equal mass 
binary neutron stars. The dashed line is the Oppenheimer-Volkoff result. Due to 
finite difference error, the numerical values systematically underestimate the mass, 
which explains why some of these values are smaller than the corresponding Oppen¬ 
heimer-Volkoff values. The insert is a blow-up of the region around the maximum 
masses. (Figure from [Baumgarte et.al. (1997|).) 


method of [Baumgarte et.al] ( |1997| , |1998b|) is based on a rescaling algorithm 
that had been used earlier for the construction of rotating stars (e.g. [I Iachisii| 
(|1986| )). Defining the point on the stellar surface that is closest to the equal 


mass companion as r A and the one that is furthest from the companion as 
r B , [Baumgarte et.al . | (|1997| , |1998b| ) specify a particular binary model by the 
maximum density and the relative separation 


z A = r A /r B . 


(9.28) 


Starting the iteration with an initial guess for the density profile, namely that 
of a spherical star, the held equations (3.27) to (3.29) can be solved. With 
a new guess for a, ^ and /3 l , the Bernoulli equation (9.24) can be evaluated 
at r A , r B and at the point of maximum density, r<> At all three points the 
density is known; at r A and r B it vanishes, and at rc it is pre-determined. 
These three relations can therefore be solved for constants D and C as well as 
the absolute separation r^^]. Given these values, the Bernoulli equation can 
be solved everywhere, yielding an updated guess for the fluid variables. The 


iteration is continued until a certain accuracy has been achieved. [Baumgarte 
et.al . | (|1997|, |1998b|) implemented this algorithm in Cartesian coordinates with 


23 In the Bernoulli equation, the gravitational fields a and i/j are rescaled in terms of 
r B , so that the latter enters the equation implicitly; see Baumgarte et.al . | (|l998b|) . 
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Mo 

Moo 

(M/R) oc 

. M 0 Disco (Jtot/M 2 ot )isco Misco 

0.0597 

0.0582 

0.05 

0.003 

1.69 

0.0578 

0.1122 

0.1066 

0.1 

0.01 

1.22 

0.1055 

0.1341 

0.1259 

0.125 

0.015 

1.12 

0.1248 

0.153 

0.1423 

0.15 

0.02 

1.05 

0.1408 

0.169 

0.1547 

0.175 

0.025 

1.00 

0.1529 

0.178 

0.1625 

0.2 

0.03 

0.97 

0.1601 


Table 9.1 

Numerical results for the orbital angular velocity D, angular momentum J and the 
ADM mass M at the ISCO, for corotating, equal mass binary neutron stars with 
polytropic index n = 1. We tabulate the individual rest mass Mq, the mass-energy 
Mqo and the compaction (M/R)^ each star would have in isolation (where R is the 
Schwarzschild radius), as well as the angular velocity MqQ, the angular momentum 
Jt 0 t/M 2 0t and the ADM mass M at the ISCO. For n = 1, the maximum rest-mass 
in isolation is M™ 3 * = 0.180. (Table adapted from Baumgarte et.al. ( 1998b[ ).) 


a full approximation storage multigrid solver. Other implementations have 
been used by [Marronetti, Mathews and Wilson] ( |1998|) ; |Uryu and Erigudn 
Gourgoulhon et.al] (|2001|) . A typical binary configuration close to the 


ISCO is shown in Figure 9.1. 


Physical units enter the problem only through the polytropic constant k in the 
polytropic equation of state (9.4). It is therefore convenient to nondimension- 
alize all equations, and present results in terms of dimensionless quantities. 
Since K n ^ 2 has units of length, the non-dimensional mass M, for example, is 
defined as 


M = k~ u/2 M, 


(9.29) 


and similar for other quantities. 


One question of interest is how the maximum allowed rest mass changes with 
separation. In Figure 9.2, the rest mass M 0 of n = 1 polytropes is plotted as 
a function of the central density p c = K n p c for different relative separations 
za■ Clearly, the maximum allowed mass increases with decreasing separation, 
suggesting that neutron stars in corotating binaries are more stable than in 
isolation. Evolutionary sequences can be constructed by keeping the rest mass 
M 0 constant. In Figure 9.2, horizontal lines therefore correspond to evolution¬ 
ary sequences. Evolving from a large separation to a smaller separation, the 
central density decreases (see also Figure 9.4 below). Similar results hold for 
other polytropic indices ( |Baumgarte et.al . | , |1998b|) . 
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The ISCO can be located by finding, for example, the minimum of the to¬ 
tal energy M along an evolutionary sequence M 0 = const. According to the 
relation 

dM = QdJ (9.30) 


(equation (33.61) in [Misner, Thorne and Wheeler| (|1973| ), compare equation 
(7.12)), a minimum in the energy should coincide with a minimum in the 
angular momentum. Values of the orbital angular velocity Af 0 Q and angular 
momentum J/M 2 at the ISCO for different values of the stellar mass M 0 are 
tabulated in Table 9.1 for n — 1 polytropes. [Baumgarte et.al .| ( |1998b|) found 
an ISCO only for sufficiently stiff equations of state, n < 1. This result is 
in qualitative agreement with various Newtonian results (see, e.g., the review 
by £ asio and Shapiro| ( |1999| )). Stars with softer equations of state are more 


centrally condensed and have more extended envelopes. Such stars come into 
contact and merge before encountering an ISCO. 


The above arguments - the increase in the maximum allowed mass and the 
decrease in the central density - suggest that neutron stars in corotating bi¬ 
naries are more stable than in isolation. A rigorous stability analysis, based 
on turning point methods, comes to the same conclusion and shows that coro¬ 
tating binaries do not encounter any instabilities until they reach the ISCO 
( |Baumgarte et.g[\ , |1998a| ). These results are in contrast to the claims of |Wilson 
and Mathew^ (|1995|) , who found in dynamical simulations that neutron stars 
collapsed to black holes individually before they reached the ISCO (“binary- 
induced collapse”; compare Section 10.2). It was later found that at least the 
size these results was erroneous and caused by a mistake in the derivation of 
the equations (|Elanagan| , |1999fe [Mathews and Wilson| . P000| ). Binary-induced 
collapse can occur for binaries made of collisionless matter (|Shapiro| , |1998| ; 
Duez et.al. 


1999|)). As seen in an inertial frame, the neutron stars in coro¬ 


tating binaries are rapidly spinning, and in fact the mass increases found in 
corotating binaries are quite similar to those found in individual stars spinning 
with the same angular velocity (see footnote [19] in [Bamiigarte et. at. \ ( |1997|) ). 
This adds one more motivation to studying irrotational binaries, which are 
more realistic and where stabilizing effects of rotation are potentially much 
smaller. 


9.3 Irrotational Binaries 


The presence of viscosity is necessary to maintain binaries in corotation. Es¬ 
timates of the viscosity in neutron stars ( |Kochanck| , |1992| [Bildsten and Cufi| 
[hT| . |1992|) show that the viscous timescale is much longer than the inspiral 


timescale, suggesting that the assumption of corotation is quite unrealistic. 
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Instead, it may be more reasonable to assume that the typical spin frequency 
of neutron stars is much less than the orbital frequency of close neutron star 
binaries. Since the emission of gravitational radiation, which drives the inspiral 
of compact binaries, conserves circulation, it is then more realistic to assume 
neutron stars to be irrotational (zero circulation), even in close binaries. 


The first relativistic formulation for irrotational binaries in quasi-equilibrium 
was given by |Bonazzola, Gourgoulhon and MarcTt| (|1997|) . Inspired by this ap¬ 
proach, less complicated formalisms were developed independently by |Asada 
(|1998|) , |Tcukolsky| ( |1998|) and [Shibata] (|1998| ) . |Gourgoulhon| ( |1998|) demon¬ 


strated that all of these formulations are equivalent. The derivation here fol¬ 
lows most closely that of |Shibata| (|1998| ). 


In general relativity, the vorticity tensor u a b is defined as 

b = P a c Pb d (Vj(ta c ) - VJJlUj)) , 


(9.31) 


where P a b = g a b + u a u b is the projection operator with respect to the fluid’s 
four-velocity u a . For irrotational binaries, the vorticity vanishes, and the quan¬ 
tity hu a can be expressed as the gradient of a potential 

hu a = V a 4>. (9.32) 


It is easy to see that with this assumption the Euler equation (9.8) is satisfied 
identically, while the continuity equation becomes 

V“ ((n//i)V a $) = 0 (9.33) 


(see, e.g., |Tcukolsky| (|1998|) ) 


Projecting these relations into E shows that for irrotational stars, the three- 
dimensional vorticity vanishes 


DjUi — DiUj = 0 , 


(9.34) 


so that 

Ui = Di$. 


(9.35) 


The Bernoulli equation (9.14) then reduces to 


h 

u l 


■Ui V i 


C, 


(9.36) 
































Velocity w.r.t corotating frame (z = 0) 


o 


o 

C\2 


a 




o 


o 

C\2 


o 



x [km] 

Fig. 9.3. The internal velocity field with respect to the co-orbiting frame in the 
orbital plane for stars of rest-mass Mq = 1.625M 0 at a coordinate separation of 41 
km, for an n = 1 polytrope with k = 1.8 x 10~ 2 Jm 3 kg -2 . The thick lines denote 
the surfaces of the stars. (Figure from |Gourgoulhon et.cd\ ( |2001| ).) 

where C is a constant. 

It is now convenient to introduce the rotational shift vector 

B i = ft + tie = r- om a . (9.37) 


From (9.10), we can express V a as 

V i = 4 -D l § - B\ 
u l h 

and the normalization u a u a = — 1 yields 
au f = (l + h^D&D^y' 2 


(9.38) 


(9.39) 


(see (9.25)). Inserting (9.38) and (9.39) into the Bernoulli equation (9.36) also 
yields the following expression for craf^J] 


au = 


-r(C + B i D& 
ah v 


(9.40) 


21 This relation can also be derived by observing that £i(hu) = £i(d<h) = d£i<h = 0 
implies that / a V a 4> = C, where C is a constant. Expressing n a in hau t = hn a u a = 
n a V a 4> in terms of l a and B l then yields (9.40); see Teukolsky ( 1998 ). 
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Fig. 9.4. The rest mass Mq as a function of central density po for separations d = 
1.3125, 1.375, 1.5 1.625, 1.75, 1.875 and 2 (thick lines from top to bottom) for 
irrotational n = 1 binary neutron stars. The dashed line is the Oppenheimer-Volkoff 
result. (Figure from Uryu, Shibata and Eriguchi ( |200C| ) .) 

The last two equations can be solved for the enthalpy h 


h 2 = a~ 2 (C + B l D&) - D i $D i $. 


(9.41) 


An equation for the velocity potential $ can be found by inserting (9.38) and 
(9.40) into the continuity equation (9.15) 


DiD'Q - B l Dj 


'C + BWj$' 


C + BWj<f> 


K 


a- 


c ^ ,c A-wb : 


a 

apo 

~h' 


(9.42) 


At the surface of the stars the density vanishes po = 0, so that regularity of 
the right hand side implies the boundary condition 

= 0. (9.43) 

surface 


'C + BWj$ 


cr 


B l - D% D iP o 


This boundary condition can also be derived by requiring that at the surface 
the fluid flow be tangent to the surface, u“V a po = 0. 

Evidently, constructing irrotational binaries is much more involved than con¬ 
structing corotational binaries. Unlike for corotational binaries, where only 
one algebraic equation (9.24) has to be solved for the enthalpy, one has to 
solve the elliptic equation (9.42) for the velocity potential <h together with the 
enthalpy equation (9.41). The boundary condition (9.43) has to be imposed 
on the surface of the star, which adds another complication since the location 
of the surface is not known a priori. 
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Fig. 9.5. Relative change of the central energy density e = po + Pi with respect 
to its value at infinite separation ej. nf as a function of the coordinate separation 
d (or the orbital frequency Q/(2n)) for corotational and irrotational constant rest 
mass sequences with Mq = 1.625 M@. The sequences were computed for a n = 1 
polytrope with k = 1.8 x 10~ 2 Jm 3 kg -2 . (Figure from Bonazzola, Gourgoulhon and 
Marckj ( |1999a|) .) 


Bonazzola, Gourgoulhon and Marck (|1999a|) and jGourgoulhon et. al\ ( [2001|) 
solved this problem with the help of a multi-domain, spectral method. In 
this method, space is covered with several patches of coordinate systems. In 
particular, the interiors of the stars are covered with spherical-type coordi¬ 
nate systems, which are constructed so that the surface of the star lies at 
a constant value of the radial coordinate. Such coordinate systems are called 
“surface-fitting” coordinates, and are very well suited for imposing the bound¬ 
ary condition (9.43). A similar algorithm, based on Newtonian simulations of 
irrotational neutron star binaries (|Uryu and Eriguchil , |1998a| Jbl), was used by 
Uryu and Eriguchil (|2000|) and |Uryu, Shibata and Eriguchil (|2000|) . [MarronettT] 
Mathews and Wilson| ( |1999|) constructed irrotational models of neutron star bi¬ 
naries in Cartesian coordinates. They simplified the boundary condition (9.43) 
by assuming that the stars are spherical, so that the gradient of the density in 
(9.43) is aligned with a radial vector. While this may be adequate as long as the 
separation between the stars are large, this approximation seems problematic 
for small separation, and may explain why the results of Marronctti, Mathew^ 
and Wilson| ( |1999|) differ from those of [Bonazzola, Gourgoulhon and Marck 
(|1999a|) ; IGourgoulhon et. at] ( |2001| ) ; |Uryu and Eriguchil (|2000|) and |Uryu, Shu 
bata and Eriguchi] ( pOOQ ) for close binaries. 


A typical binary configuration and its internal velocity field is shown in Fig¬ 
ure 9.3. The maximum allowed mass of neutron stars in irrotational binaries 
can be found by hireling the mass as a function of central density for fixed 
separation. Results of |Uryu, Shibata and Eriguchil (|2000|) are shown in Fig¬ 
ure 9.4. This demonstrates that, as in corotational binaries, the maximum 














































































Mo 

Moo 

(M/R) oo 

MoOcusp 

(J to t/M t 2 ot )cusp 

Mcusp 

0.112 

0.107 

0.1 

0.0106 

1.09 

0.105 

0.130 

0.122 

0.12 

0.0144 

1.02 

0.121 

0.146 

0.136 

0.14 

0.0187 

0.971 

0.135 

0.166 

0.153 

0.17 

0.0263 

0.919 

0.150 

0.175 

0.160 

0.19 

0.0320 

0.895 

0.157 


Table 9.2 

Numerical results for the orbital angular velocity f2, the angular momentum J and 
the ADM mass M at cusp formation, for irrotational binary neutron stars with 
polytropic index n = 1. We tabulate the individual rest mass Mq, the mass-energy 
Mqq, the compaction ( M/R) 00 each star would have in isolation, as well as the 
angular velocity MqQ, the angular momentum J to t/M 2 ot and ADM mass M at cusp 
formation. For n = 1, the maximum rest-mass in isolation is M™ ax = 0.180. (Table 
adapted from|Uryu, Shibata and Eriguchj (2000|).) 


mass increases with decreasing separation. However, comparing with Figure 
9.2, we find that the increase in maximum mass is smaller for irrotational 
binaries than for corotational binaries (especially taking into account that the 
fairly coarse resolution results of [Baumgarte et.al . | (|1998b|) underestimate the 
masses in Figure 9.2). This result is not surprising, since neutron stars in 
corotational binaries are rotating with respect to the rest frame of the binary, 
which by itself increases their maximum mass (e.g. |Cook, Shapiro and TcukoT| 
sky| (|1994|) ). It is also evident from Figures 9.4 and 9.2 that while the density 
along evolutionary sequences Mq = const of irrotational binaries decreases 
with decreasing separation, the decrease is less than for corotational binaries. 
This can be seen more clearly in Figure 9.5, which demonstrates that while 
it is possible that the central density increases with decreasing separation by 
a very small amount for intermediate separations in irrotational binaries, it 
certainly decreases for small separations where tidal deformation dominates 
over any other effects. Similar results were found by |Uryu and Eriguchi| ( |2000|) 
and |Uryu, Shibata and Eriguchi| (|2000|) . While these findings have no imme¬ 
diate implications for the stability of neutron stars in irrotational binaries, 
they offer no evidence for an instability as initially reported by [Wilson and 
Mathews) flT995l) . 


While evolutionary sequences of corotational binaries end either at the ISCO 
or at contact, irrotational sequences typically end with the formation of a cusp 
before they reach the ISCO ([Bonazzola, Gourgoulhon and Marck[ , |1999a| ; [Uryvi 
and Eriguchi). |200(J|; |Uryu, Shibata and Hrigucln|, j2000|). As analyzed by |Uryu 


Shibata and Eriguchi| (|2000|) , this cusp corresponds to an inner Lagrange point, 


across which neutron stars will transfer mass. Such configurations are likely to 


form dumbbell-like structures, with mass flowing between the two stars. Uryu. 


Shibata and Eriguchi| (|200q) find that only binaries with very stiff equations 


90 











































































of state (n < 2/3) reach an ISCO before they form a cusp, while binaries with 
softer equations of state form a cusp first. Numerical values for irrotational 
n — 1 binaries at cusp formation are listed in Table 9.2. Comparing with 
the ISCO values of corotational binaries in Table 9.1, one finds that the cusp 
and ISCO occur at very similar frequencies. The corotational binaries have 
more angular momentum, because the individual stars carry a spin angular 
momentum in addition to the orbital angular momentum of the binary. We 
also find that the binding energy (M — M oo )/M 0 of corotational binaries is 
slightly larger than for irrotational binaries, because the ADM mass M of the 
former include the additional spin kinetic energy of the individual stars (see 
also the discussions in |I)uez et7aT\ ( |2002| ) and Section 10.4). More numerical 
results and a detailed discussion of the ISCO and cusp formation in irrotational 
binaries can be found in |Uryu, Shibata and Eriguchi| (|2000|) . 


10 Dynamical Simulations of Binary Neutron Stars 


In this Section we present approaches and results for dynamical simulations of 
binary neutron stars. We discuss formulations of relativistic hydrodynamics in 
Section 10.1, the conformal flatness approximation in Section 10.2, fully self- 
consistent simulations in Section 10.3, and a quasi-adiabatic approximation of 
the slow inspiral just outside of the ISCO in Section 10.4. 


10.1 Relativistic Hydrodynamics 


Since numerical hydrodynamics in general relativity has been the subject of a 
recent review (|Font| , |2000| ) we will provide only a very brief overview. 

The equations of relativistic hydrodynamics can be derived from the local 
conservation laws of the stress-energy tensor T ab (9.7) 

\7 a T ab = 0 (10.1) 


and the continuity equation (9.9) 

V a (p 0 w a ) = 0, 


( 10 . 2 ) 


where p 0 is the rest mass density and u a the fluid four-velocity. For a perfect 
fluid, the stress-energy tensor is (9.2) 

T ab = (p 0 + Pi + P)u a u b + Pg ab (10.3) 


91 


















where pi is the internal energy density and P the pressure. 


The equations of motion can be derived from (10.1) and (10.2) in various 
different ways, depending on how the basic dynamical fluid variables are de¬ 
fined. One might be tempted to adopt a projection of T ab with respect to the 
normal vector n a , which yields the matter variables p (2.26) and j l (2.29) as 
they appear in the ADM equations. It turns out, however, that projections 
with respect to the fluid four-velocity u a yields equations that are simpler and 
closer to the Newtonian equations of hydrodynamics. 


Wilson ( |1972| , see also |Hawlcy, Smarr and Wilson! (|1984a| .[b|); [Wilson and Ma.th- 
ews| (|1989|) ; |Wilson, Mathews and Marronetti] ( |1996|) ), who pioneered the held 
of relativistic hydrodynamics in multi-dimensions, used the density definition 
D = WPq, where W is the Lorentz factor between the fluid four-velocity u a 
and normal observers n a 


W = n a u a = era*. 


(10.4) 


The equation of continuity (10.2) then becomes 


d t (y/jD)+d j (y/jDv^)= 0 . 


(10.5) 


Here v l = u l /u l is the transport velocity with respect to a coordinate observer, 
and 7 is the determinant of the spatial metric 7 ^. Defining the momentum as 
S a = pohWu a , the spatial components of (10.1) reduce to the Euler equation 


1 


dtiyfiSi) + 


1 

7^ 


3;(VW) 


-adiP 


SgS b 
2 S° 


dig 


ab 


( 10 . 6 ) 


The time component S° is found from the spatial components through the 
normalization u a u a = —1. With the energy density E = Wpi, the timelike 
component of ( 10 . 1 ) yields the energy equation 

dt(VlE) + ^(vW) = ~P (dt{y/W) + diiy/Wv*)) ■ (10.7) 




] 


Shibata (p.999a| , see also |Shibata, Baumgarte and Shapiro| (|1998| ) ) modihed 
this scheme by absorbing the determinant into the definition of the density 
p * = -\fyWpo, for which the equation of continuity becomes 


d t p * + djip-^v 1 ) = 0 . 


( 10 . 8 ) 


Also, for gamma-law equations of state 

p = (r -1 )pi 
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(10.9) 








































the source term on the right hand side of the energy equation (10.7) can be 
eliminated 


<9 t e* + dj(e*v j ) = 0 ( 10 . 10 ) 

i /r 

if the energy variable is defined as e* = sJfiWpfi . Defining Ui = hu the 
Euler equation (10.6) becomes 


d t {p*Ui) + dj(p*UiV J ) 


—ae e<l> diP — p* (Whdia — Ujdift 
ae-^UjUk jk 2 ah(W 2 - 1) 
+ 2 Wh a W 



( 10 . 11 ) 


where we have expressed the right hand side of ( 10 . 6 ) in terms of the metric 
quantities introduced in the BSSN formalism (Section 4.3; in particular = 
e 6 0 )_ qq ie transport velocity v l can be found from 


-/? + 


CK7 lj Uj 

Whe 4 P 


( 10 . 12 ) 


and W can be found from the normalization (9.25), which can be expressed 
as 


W 2 


Y j UjUj 

e^ 



re 


-2 


p*(We^) 


r-i 


(10.13) 


Example 10.1 For static configurations we have Ui = 0 and hence W = 1 
and, in the BSSN formulation, p* = e^po. The Euler equation (10.11) then 
reduces to 


adiP + pohdia = 0 , 


(10.14) 


which can be combined with (9.6) to yield equation (9.19) of Example 9.1 for 
hydrostatic equilibrium. 

Other groups have cast the equations of relativistic hydrodynamics into the 
flux-conservative form 

d t U + diP = S (10.15) 


where U is the state vector containing the so-called primitive fluid variables, 
T is the flux vector, and the source vector S does not contain any derivatives 
of the fluid variables (see, e.g., [Marti, Ibanez and Mirallcsj ( [L991|) , [Font] ( j2000| ) 
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and references therein). [Font et.aT . ] (|2000| , [200 1|) have implemented such a 
scheme with 


U 


1 

% 



y/Wpo 

blPohW 2 Vj 

K ^f(p 0 hW 2 -P-Wp 0 ) / 


(10.16) 


where v l = u l /W + (3 1 /a is the velocity of the fluid with respect to a normal 
observer, and where 7 is the determinant of the spatial metric 7 ^. The flux 
vector T is then given by 


F 


(a# - :f)b 

(oiv 1 - /3 l )Sj + a^PyP5 l j , 
^ ( av l — /3 l )f + oiyjFPP j 


(10.17) 


and the source vector S by 

( n \ 


5 = 


0 


a^T ab g hc ^Y 


[a^(T a0 d a a-aT ab ^r° ab ) 


(10.18) 


The virtue of these schemes is that the local characteristic structure can be 
determined, which is crucial for the implementation of high-resolution shock- 
capturing schemes (see below). 

Once the equations have been brought into a particular form, a numerical 
strategy for their numerical implementation has to be chosen (see [Font| ( p000|) ). 
This is true for all equations in this article, but hydrodynamics poses the addi¬ 
tional challenge that shocks may appear. In a shock, hydrodynamic quantities 
develop discontinuities when perfect fluids are assumed, and macroscopic fluid 
flow is converted into internal energy. Neither one of these phenomena can be 
captured by traditional finite difference schemes, which therefore have to be 
modified. 


Von [Neumann and Richtmyer| (|1950| ) suggested an artificial viscosity scheme, 
in which an artificial term Q is added to the pressure P. They showed that 
the jump conditions across shocks are well satisfied for a Newtonian, one- 
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dimensional fluid flow if Q is defined as 


( — p(kAx) 2 (d x v) 2 if d x v < 0 

Q= \ (10.19) 

I 0 otherwise, 


where Ax is the grid resolution and k an adjustable constant of order unity. 
The effect of the artificial viscosity is to spread out the shock over approxi¬ 
mately k gridzones. It appears in the source terms of the Euler equation (10.6) 
and the energy equation (10.7), and, as desired, converts the kinetic energy 
of bulk fluid flow into internal energy in accord with the Rankine-Hugoniot 
jump conditions. 


Artificial viscosity schemes have been generalized for applications in general 
relativistic hydrodynamics by numerous authors (including [May and Whitc| 
(|1966|) ; |Wilson| (|1972|) ; |van Riper] (|1979|) ; [Shapiro and Tcukolskyj (|1980|) ; |Baum- 
gartc, Shapiro &; Tcukolsky] ( |1995|) ; [Shibata| (|1999a|) ). These schemes have the 
virtue of being quite robust and very easy to implement. However, it has been 
shown that artificial viscosity schemes do not work well for ultra-relativistic 
fluid flow (|N orman and Winkler], |1986|). 


Alternatively, high-resolution shock-capturing (HRSC) schemes can be em¬ 
ployed. These schemes are based on the idea of treating all fluid variables 
as constant in each grid cell. The discontinuity at each grid interface can then 
be viewed as the initial data for a Riemann shock tube problem, which can be 
solved either exactly or approximately. These schemes rely on the local char¬ 
acteristic structure of the equations for the solution of the Riemann problem, 
and are therefore used with equations in flux-conservative form, ( Font et7aT\ , 
2000 | , |200 1| ) . By construction HRSC are capable of handling shocks and do not 
need additional artificial viscosity terms. In various applications, HRSC have 
been found to be more accurate than artificial viscosity schemes, in particular 
for highly relativistic flows. For more information on HRSC schemes in rela¬ 


tivistic hydrodynamics see the review articles by Marti and Muller (1999) and 
Font] (goog). 


It is also possible to avoid finite difference methods altogether, and to solve 
the equations of relativistic hydrodynamics with either smoothed particle hy¬ 
drodynamics (SPH) or spectral methods. SPH methods are quite common in 
Newtonian hydrodynamics, and have also been applied in relativistic hydro¬ 
dynamics in stationary spacetimes (see |Font| ( |2000|) for references). Quite re¬ 


cently, SPH methods have also been used for the modeling of binary neutron 
star coalescence in post-Newtonian theory ( [Faber and Rasic| , |2002|) as well as 
in full general relativity ( |Occhslin, Rosswog and Thiclcmann| , |2002|) . Spectral 
methods (see [Bonazzola, Gourgoulhon and Marck| (|1999b|) for a review) have 
been used quite commonly for the solution of the elliptic equations in ini- 
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tial value problems (e.g. [Boiiazzola, Gourgoulhon and Marck| (|1999a|) ; |Grand- 


clement, Gourgoulhon and Bonazzolal ( |2002| ) ; IGourgoulhon et.al\ ( |2001|) ) and 


recently also for the dynamical evolution of gravitational (vacuum) fields ( |lvid-| 
dcr, School and Tcukolskyl , [200 1|) . For hydrodynamic simulations, in particular 


in three dimensions, spectral methods are less common because of difficulties 
in treating the discontinuities in shocks, which lead to spurious oscillations 
(Gibbs phenomenon). 


10.2 The Wilson-Mathews approximation 


Wilson and Mathews ([1989], |1995| , also |Wilson, Mathews and Marronetti| ( |1996| ) ) 
pioneered the conformal flatness approximation to the simulation of neutron 
star binaries. In this approach, one assumes that the dynamical degrees of 
freedom of the gravitational fields, i.e. the gravitational radiation, play a neg¬ 
ligible role for the structure of binary neutron stars. To simplify the metric 
and the held equations, [Wilson and Mathews) (|1989[ |1995|) therefore suggest 
assuming that the spatial metric is conformally hat, so that the spacetime 
metric takes the simplified form 


ds 2 = —a 2 dt 2 + ip 4 rjij(dx l + (5 l dt){dx^ + ftdt). 


( 10 . 20 ) 


Wilson and Mathews (|1989[ |1995|) further assume that the spatial metric re¬ 
mains conformally hat at all times. The traceless part of equation (2.46) then 
has to vanish, which yields 

A ij = ^-(L/3)W (10.21) 


Here kW is the traceless part of the extrinsic curvature, and the vector gradient 
L is defined in (3.10). 

Equation (10.21) can be inserted into the momentum constraint (2.45), re¬ 
sulting in an equation for the shift [3 l . The conformal factor ijj can be found 
from the Hamiltonian constraint (2.44), and to determine the lapse a |Wilson| 
and Mathcwsl (|1995|) adopt maximal slicing (Section 5.2). If one also adopts 
the conformal rescaling (3.7), these equations reduce to the thin-sandwich 
equations (3.27) through (3.29), which now completely determine the metric 
( 10 . 20 ). 

All unknowns in the metric (10.20) are determined by elliptic equations, and 
in this sense dynamical degrees of freedom have been removed from the grav¬ 
itational fields. In this approach, one solves an initial value problem at each 
instant of time, as opposed to dynamically evolving the gravitational fields. 
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While one may worry about the accuracy of this approximation (see also Ap¬ 
pendix C), it greatly simplifies the held equations and allowed [Wilson and| 
Mathews| (|1995|) to perform the first fairly detailed simulations of binary neu¬ 
tron stars. In this approach, the time step is limited by the matter sound speed 
and not the light speed, so it can be much larger than in fully self-consistent 
algorithms. 


In these simulations, |Wilson and Mathews 


combined equations (3.27) 
through (3.29) for the metric (10.20) with Wilson’s formulation of the rela¬ 
tivistic equations of hydrodynamics (10.5) through (10.7). In each timestep, 
one first evolves the matter variables, and then solves the held equations for 
the metric (10.20) with the new density distribution as sources. If desired, the 
gravitational wave emission can be estimated with the quadrupole formula. 


The initial results of Wilson and Mathews| (|1995|) showed a so-called crushing 
effect, in which the neutron stars were compressed to very high densities even 
far outside of the ISCO, and ultimately collapsed to two individual black holes 
(“binary-induced collapse”). Since this effect contradicted expectations from 
Newtonian and post-Newtonian theory as well as relativistic quasi-equilibrium 
results (which predict that the density should be reduced by tidal deforma¬ 
tions, see Sections 9.2 and 9.3), their findings spurned a fairly intense debate 
until it was discovered that at the very least the size of the effect was caused by 
an inconsistency in the derivation of the equations (|Flanagan|, |1999| ; [Mat he w^ 
and Wilson| , |2000|) . 


The unfortunate consequence of this debate is that the assumption of confor¬ 
mal flatness itself was often incorrecty blamed for the spurious results. While 
this approximation is obviously not suitable for every situation, it may be very 
good for the modeling of compact objects in some regimes, where its errors are 
probably much smaller than other sources of error (see also the discussion in 
Appendix C). It has recently been adopted by Pcchslin, llosswog and Thiclc- 
man 3 (Eli. who combined it with an SPH method to solve the equations of 
relativistic hydrodynamics. 


10.3 Fully Self-consistent Simulations 


Several groups have launched efforts to self-consistently solve the equations 
of relativistic hydrodynamics together with Einstein’s equations, and model 
the coalescence and merger of binary neutron stars ( |Uohara and Nakamura . 
1997| ; [Baumgarte et.al . | . 1999| ; [Font et.aJ\ , |2000[ [200 1| ). The first successful 


simulations of binary neutron star mergers are those of [jhibata and Uryh 
(|20004 see a l so |Shibata| (|1999a|) ; |Shibata and Uryh| (|2002 )), which we will 
focus on in this Section. 
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Model 

Pmax 

Mo 

M 

J/M 2 

remnant 

11 

0.0726 

0.261 

0.242 

0.98 

neutron star 

12 

0.120 

0.294 

0.270 

0.93 

black hole 

13 

0.178 

0.332 

0.301 

0.88 

black hole 


Table 10.1 

Summary of the initial data for the coalescence simulations of n = 1 irrotational 
binary neutron stars by Shibata and Uryd| (|2000a ). Here Mq, M and J are the total 
rest mass, mass and angular momentum of the binary. In these dimensionless units, 
the maximum allowed rest mass of an isolated, non-rotating star is M“ ax = 0.180. 
(Table adapted from [Shibata and Uryu| (|2000a|) .) 


Shibata and Uryu ( p 000 a|) adopt a gamma-law equation of state (10.9) with 
index T = 2 [n = 1). They solve the equations of relativistic hydrodynamics 
in the form (10.8) to (10.11), and Einstein’s equations in the original form 
( |Shibata and Nakamura| , |1995|) of what is now often referred to as the BSSN 
equations (Section 4.3). They use “approximate maximal slicing” (Section 5.2) 
to specify the lapse a and “approximate minimal distortion” (Section 5.4) to 
determine the shift (3 l . They also add a radial component to the shift vector 
to avoid grid-stretching in collapse situations. 


As initial data, |Shibata and Uryh| ( [2000a|) prepare equal-mass polytropic 
(n = 1 ) models of binary neutron stars in quasi-equilibrium with both coro- 
tational (see Section 9.2) and irrotational (Section 9.3) velocity profiles. For 
both velocity profiles they prepare three different models with individual stel¬ 
lar masses ranging from about 75% to 100% of the maximum allowed mass 
of non-rotating stars in isolation. For corotating models, |5hibata and Uryh 
(|2000a|) adopt the contact models (za = 0 in the parametrization (9.28)) as 
initial data, which are fairly close to the ISCO (see Section 9.2). Irrotational 
sequences terminate at cusp formation (see Section 9.3), which is still outside 


of the ISCO. [Shibata and Uryu| 


therefore adopt the cusp model as ini¬ 


tial data, and induce collapse by artificially reducing the angular momentum 
by about 2.5%. Since irrotational velocity profiles are probably more realis¬ 
tic, we will discuss their models (II), (12) and (13) for irrotational binaries. 
The initial data for these three models are summarized in Table 10.1 (more 
information can be found in Table 1 of |Shibata and Uryii| ( [2000a|) ). 


The Erst simulation, model (II) in IShibata and Uryu| ( |2000a| ), is for a binary 
of total rest mass M 0 = 0.261 Q The angular momentum of the initial data 
is J/M 2 = 0.98, where M is the total gravitational mass, and is hence smaller 
than the Kerr limit J/M 2 = 1. We show snapshots of density contours in 
Figure 10.1. 


25 We have converted the units of Shibata and Uryu ( [2000a ) into the same dimen¬ 
sionless units adopted in Section 9 for easier comparison. 


98 














































t= 1.41P 


t= 1.81P 




X/ M g 


X/ M g 


Fig. 10.1. Snapshots of density contours for p* (see equation (10.8)) and the ve¬ 
locity field (v x ,v y ) in the equatorial plane for the coalescence of an irrotational 
binary of total rest mass Mo = 0.261 (Model (II) in Shibata and Uryu ( p000a|) ). 
Time is measured in terms of the orbital period P. The contour lines denote den¬ 
sities p* / p * max — 10 0 3j with p* max = 0 .255 and j = 0.1, 2, • • • , 10. (Figure from 
Shibata and Uryij ( 2000a| ).) 


In contrast to the coalescence of corotational binaries, no significant spiral 
arms form during merger, and hardly any matter is ejected. It is quite sur¬ 
prising then that the remnant settles down to a near-equilibrium neutron star 
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X / M g , Y / M g 


Fig. 10.2. The angular velocity along the x-axis (solid line) and the y-axis (dotted 
line) at t = 1.81P or b for model (II). Here M g is the gravitational mass M. (Figure 
from [Shibata and Uryu| ( |2000a| ).) 


and does not collapse to a black holcp^], even though its rest mass exceeds the 
maximum allowed rest mass of a spherical, non-rotating star by about 45%. 
Shibata and Uryil ( ^UUUa ) find only small amounts of shock heating, which 
rules out thermal pressure as providing the extra support. As we have pointed 
out, the angular momentum J/M 2 is smaller than the Kerr limit and there¬ 
fore cannot prevent black hole formation. Uniform rotation can increase the 
maximum allowed rest mass for T = 2 polytropes by only about 20% (|Cook 


Shapiro and Tcukolsky| , |1994|) . However, |5hibata and Uryh| (|2000a|) find that 
the core is differentially rotating as opposed to uniformly rotating, as shown 
in Figure 10.2. 

It is quite intuitive that differential rotation may significantly increase the 
maximum mass ([Baumgartc, Shapiro fc Shibata| , |2000|) . For uniform rotation, 
the angular velocity, and hence the centrifugal force which balances the grav¬ 
itational force to increase the maximum mass, is limited by the Kepler limit 
at the equator, above which matter there would no longer be gravitationally 
bound (the “mass-shedding limit”). For differential rotation, the core may ro¬ 
tate faster than the equator, and may further increase the maximum mass 
without violating the Kepler limit at the equator. Paumgarte, Shapiro fc SET 


bata| ( pOOOl ) have found that differential rotation is very effective in increasing 
the maximum mass, and have found mass increases of over 60% even for very 
modest degrees of differential rotation. 


The remnant formed in this simulation is supported by differential rotation, as 


26 At least not on a dynamical timescale - see below. 
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X/ M g 


X/ M g 


Fig. 10.3. Same as Figure 10.1, but for a binary of total rest mass Mq = 0.332 
(Model (13) in [Shibata and UryO ( 2000a )). Contour lines denote densities 
p*/p* max = 10 -0 ' 3j with p* max = 0.866 and j = 0,1, 2, • • • , 10. The dashed line 
in the last snapshot is the circle r = 3 M which encloses over 99% of the total rest 
mass. The thick solid line at r ~ M denotes the location of the apparent horizon. 
(Figure from Shibata and Uryu (|2000a| ).) 


the thermal pressure support is minimal while its rest mass exceeds the maxi¬ 
mum allowed rest mass by about 45%. Similar results were found more recently 
by |Oechslin, Rosswog and Thielemann| ( |2002| ) , who used a completely indepen¬ 
dent numerical method (namely an SPH method to model the hydrodynamics 
and the Wilson-Mathews approximation to model the gravitational fields). 
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Fig. 10.4. Gravitational wave amplitudes h + and h x as a function of retarded time 
for the irrotational binary neutron star models (E) and (H) of [Shibata and UryQ 
(2002). Model (E) corresponds to a slightly smaller mass then (II) and results in 
a differentially rotating neutron star, while model (H) is similar to model (12) and 
results in a black hole. (Figures from Shibata and Uryii| (|2002[).) 


This result is quite surprising, since up to then it was usually assumed that 
such massive neutron stars collapse to black holes on a dynamical timescale, 
(although earlier Newtonian merger calculations already foreshadowed this 
result; see |Rasio and Shapiro| ( |1999| ) and references therein). Instead, it is pos¬ 
sible that such “hyper-massive” remnants are dynamically stable and collapse 
to black holes on a secular timescale, after some dissipative mechanism, for 
example viscosity or, more likely, magnetic fields, has brought the star into 
more uniform rotation (|Baumgartc, Shapiro fc Shibata| , [20001 ; |Shapiro| , |2000|) . 


Both models (12) and (13), for which the total rest mass exceeds the maximum 
allowed mass by 63% and 85% (see Table 10.1), form a black hole within a 
dynamical timescale upon merger. In Figure 10.3 we show snapshots of the 
density contours and velocity profiles for the more massive model (13). In 
the last frame of Figure 10.3 the thick solid line denotes the location of an 
apparent horizon (see Section 6.1), indicating the formation of a black hole. 

I I __, 

Shibata and Uryu (|2000a|) find fairly similar results for corotational binary 
models. Probably the most significant difference is that corotational binaries 
have more angular momentum in the outer parts of the binary, which leads to 
the formation of spiral arms during the coalescence. The spiral arms contain 
a few percent of the total mass, and may ultimately form a disk around the 
central object. 


In the simulations of |Shibata and Uryh| (|2000a|) one of the largest limita¬ 
tions on the accuracy were the outer boundaries, which, because of limited 
computational resources, had to be imposed well within a wavelength of the 
gravitational radiation from the binary (at about Agw/3). This means that the 
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waves were extracted without being in the radiation zone, which necessarily in¬ 
troduces error. After having gained access to a more powerful supercomputer, 
Shibata and Uryu| ( j2002| ) therefore repeated these calculations on computa¬ 
tional grids that extend to about a gravitational wavelength. Qualitatively, 
these improved results are very similar to their earlier ones, although the on¬ 
set of black hole formation shifts to slightly smaller masses. Most strongly 
affected by these improvements are the gravitational wave forms. In Fig 10.4, 
we therefore show examples from these improved simulations for models that 
are similar to model (II), leading to a differentially rotating neutron star, and 
model (12), leading to a black hole. 


While the simulations of [Shibata and IJrvul (|2000a|, |2002|) are pioneering, it 
would be desirable to confirm their findings with independent simulations with 
fully self-consistent codes. Many aspects of the simulations could also be im¬ 
proved in the future, including the setup of the initial data (eliminating the 
need for an artificial reduction of the angular momentum), the extraction of 
gravitational radiation, and the handling of the hydrodynamics. For exam¬ 
ple, coalescing irrotational neutron stars form a vortex sheet at the contact 
surface that is Kelvin-Helmholtz unstable on all wavelength (see, e.g., [Rasio 
and Shapiro| (|1999|) and references therein). Reliably simulating such sheets is 


quite challenging, and is likely to require more sophisticated algorithms than 
artificial viscosity schemes (compare Section 10.1, [Font et.aT . | (|2000[ [200 1|) ). 
Lastly it should be emphasized that these simulations assume polytropic stars 
governed by a gamma-law equation of state, which is an idealization. In re¬ 
ality, a number of factors, including effects of a realistic nuclear equation of 
state, magnetic fields and neutrino transport, may play an important role in 
the coalescence of binary neutron stars. It will probably be a while until all 
these can be incorporated in fully dynamical and self-consistent simulations. 


10.4 The Quasi-Adiabatic Inspiral of Binary Neutron Stars 


As we have discussed in Section 1, it is possible that other means of model¬ 
ing binary neutron stars, in particular PN point-mass techniques, break down 
somewhat outside of the ISCO, when finite-size and relativistic effects become 
important. It is hard to imagine that fully hydrodynamical numerical calcu¬ 
lations would be able to follow the inspiral reliably from such a point outside 
of the ISCO through many orbital periods to the onset of instability at the 
ISCO, followed by plunge and merger. Such calculations would accumulate sig¬ 
nificant amounts of numerical error and would be computational prohibitive. 
This leaves a gap between the regimes that PN and fully numerical calcula¬ 
tions can model. Filling this gap for the late inspiral, immediately prior to 
plunge, therefore requires an alternative, approximate approach (in the case 
of binary black holes, this problem has been called the “intermediate binary 
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Fig. 10.5. The final hundreds of cycles of the inspiral waveform h+ or h x on the 
axis of rotation as a function of retarded time and cycle number for an irrotational 
binary neutron star system. Also indicated is the gravitational strain h for a binary 


of total rest mass Mq = 2 x 1.5 Mq at a separation of 100 Mpc. (Figure from Duez 

B') 


et. al. 


black hole” problem ( Brady, Creighton and Thorne , 19981 )). 


Several different such approaches have been suggested ( [Blackburn and 1~ 


tweiler 

, 1992; 

Brady, Creighton and Thorne|, |1998|; |Laguna, 

1999 

; W helan 

and Romano, 

1999|; Whelan, Krivan and Price|, 

2000 

; Duez, Baumgarte and 


Shapiro| , |2001| ; |Duez et.~aJ\ |2002| ; |Shibata and Uryu| , |2001|) . Here we will fo¬ 
cus on the “hydro-without-hydro” approach (see also [Baumgartc, Hughes and 
Shapiro! ( |1999| ) ) adopted by puez, Banmgarte and Shapiroj (|2001| ) and puez 
et.al. (|2002| ) (see also |Yo, Baumgarte and Shapiro| (|2001a|) , who illustrated 


and calibrated this approach in a scalar gravitation model problem). 


This approach approximates the binary orbit outside of the ISCO as circular, 
and treats the orbital decay as a small correction. For each binary separa- 
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Fig. 10.6. A match of the late quasi-equilibrium inspiral wavetrain and the plunge 
and merger waveform for a binary of total rest mass Mq = 1.62M™ ax . The merger 
of these binary results in immediate black hole formation (compare Section 10.3). 
Here h x = r s /Mo^y xy , where r s is the distance to the source. For Mo = 2 x 1.5 Mq 
and a distance to the source of 100 Mpc, the maximum amplitude of the metric 
perturbation is ~ 5 x 10~ 21 . (Figure from Duez et.al. ( 2002 ).) 


tion, the matter distribution can then be determined independently by the 
quasi-equilibrium methods of Section 9. These matter profiles, which satisfy 
the equations of quasi-stationary equilibrium, can then be inserted as matter 
sources into Einstein’s equations. Evolving the gravitational fields yields the 
gravitational wave signal and luminosity, and hence the rate at which the sys¬ 
tem loses energy at each separation. Combining this luminosity dM/dt with 
the derivative of the binding energy M with respect to separation r yields the 
inspiral rate 


dr dM / dt 
dt dM / dr 


( 10 . 22 ) 


Integrating this equation yields the separation as a function of time, and ac¬ 
cordingly the entire gravitational wave train. 


Duez, Baumgarte and Shapiro ( |2001| ) implemented such a scheme for coro- 
tational binaries, based on the quasi-equilibrium models of [Baumgarte et.al.\ 
(|1998b| ) (see Section 9.2). In |Duez et.al] (|2002|) , these results were compared 
with those for an irrotational sequence, based on the models of Uryu, Sim 
|bata and Eriguchl| ( |2000| ) (see Section 9.3). These simulations are illustrative 
only due to the small outer boundary radius which necessitated gravitational 
read-off inside the wave zone (compare Section 10.3). However, the method is 
quite promising. In Figure 10.5 we show such a gravitational wavetrain for an 
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irrotational binary. 


For a given separation, the gravitational wave luminosity dM/dt is very similar 
for corotational and irrotational models[ ?7 ~|. but since the binding energy of 
corotational models includes the spin kinetic energy of the individual stars, 
the absolute value of the corotational binding energy and hence \dM/dr\ is 
smaller than that of irrotational binaries. According to (10.22) this makes 
\dr/dt\ larger, meaning that the inspiral of corotational binaries proceeds faster 
than that of irrotational binaries. 


Duez et.al. ( |2002| ) also point out that the entire gravitational wavetrain, from 
the slow inspiral to the ISCO and the subsequent plunge and merger, can be 
constructed by matching results from an quasi-adiabatic approximation of the 
inspiral and a dynamical simulation of the coalescence. An example, showing 
the last 13 orbits of the inspiral together with the plunge and merger leading 
to black hole formation (as obtained by [Shibata and Uryh| (|2000a|) ), is shown 
in Figure 10.6. 


11 Summary and Outlook 


A search of the recent literature reveals the impressive progress that numerical 
relativity and the modeling of compact binaries has made in the past few years. 
New formulations of the initial value problem (Section 3) and the evolution 
equations (Section 4), new coordinate conditions and their implementations 
(Section 5) as well as new diagnostics (including the horizon finders described 
in Section 6) have led to several advances in the simulation of black holes 
and neutron stars. Recent breakthroughs include initial data for binary black 
holes (Section 7) and neutron stars (Section 9), evolution calculations for single 
and binary black holes (Section 8) and simulations of inspiraling and merging 
binary neutron stars (Section 10). Some results of astrophysical interest are 
density profiles of relativistic binary polytropes, the location of the ISCO 
in binary black holes and neutron stars, the stability properties of neutron 
stars in close binaries, cusp formation in irrotational binary neutron stars, the 
stabilizing effect of differential rotation in remnants of binary neutron star 
coalescence, and the first preliminary gravitational waveforms from coalescing 
binary neutron stars and black holes. 

One of the long-standing goals of numerical relativity is the simulation of the 
coalescence of binary black holes. While this goal has not been achieved yet, 

27 As one might expect, since the gravitational wave emission is dominated by the 
matter density, which is fairly similar for the corotational and irrotational binaries, 
while matter current distributions play a less important role. 
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some preliminary calculations have been performed (Section 8.2). Moreover, 
with recent progress in the evolution of single black holes (Section 8.1), the 
simulation of binary black holes seems much more feasible than it did only a 
few years ago. 


It is likely that most of the numerical results reported in this article will be 
revisited and improved. Some of these possible improvements have already 
been discussed in the corresponding Sections. For binary black hole initial 
data (Section 7) it would be desirable to gain a clearer understanding of why 
different approaches lead to different values of the ISCO. Specific calcula¬ 
tions, some of which have already been initiated, include improvements of the 
thin-sandwich approach of |Grandclement, Gourgoulhon and Bonazzola ( |2002|) 
and the construction of binaries in circular orbit using Kerr-Schild or post- 
Newtonian background data, which would avoid the assumption of spatial 
conformal flatness. Avoiding this assumption may also be interesting for bi¬ 
nary neutron star initial data (Section 9), in addition to extending sequences 
of irrotational binaries beyond cusp formation (see Section 9.3). The latter 
would provide better initial data for dynamical simulations of binary neutron 
star coalescence (Section 10.3), which could also be improved by using more 
sophisticated hydrodynamics algorithms. Other future improvements will be 
the inclusion of more realistic nuclear equations of state, magnetic fields, and 
possibly neutrino transport. Most of the results that we have discussed adopt 
the assumption of equal mass binaries, which eventually will also have to be 
generalized 28 . 


Almost all numerical results discussed in this article would benefit from in¬ 
creased accuracy, and most will need to be significantly more accurate before 
they can be useful in gravitational wave catalogs as discussed in Section 1. 
On uniform grids, two competing sources of numerical error are finite differ¬ 
ence errors (due to coarse numerical grid resolution) and the location of outer 
boundaries (if they are imposed at too small a separation). Given the avail¬ 
able computational resources and hence the size of the computational grid that 
can be afforded, a more or less suitable compromise has to be chosen between 
grid resolution and distance to the outer boundaries. In current simulations, 
this compromise still leads to significant error, in particular for gravitational 
wave forms (see Sections 10.3 and 10.4). Two improvements may be possible: 
instead of uniform grids, nested grids or adaptive mesh refinement (AMR) 
may be employed. While such techniques have been used in many other fields 
of computational physics as well as in lower-dimensional numerical relativity 
(see, e.g. |Choptuik| ( |1993| )), they have been used sparingly in 3 + 1 numeri¬ 
cal relativity so far (e.g. jErugmarm] (|1996| , |1999|) ; [New et.al\ ( |2000| ); |Jansen 


^_Note, however, that the mass of most neutron stars in neutron star binaries is 
close to 1.4 Mq ( [Thorsett and Chakrabarty , 1999), so that the assumption of equal 
mass seems quite adequate. 
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et.al. (|20Q2|) ). In addition, the handling of the outer boundaries and gravita¬ 
tional radiation extraction could be refined. Some recent efforts include Bishopl 


et.al . | ( |1996| ) ; |Abrahams et.al .| (|1998| ) ; |Fricdrich and N agy| (|1999|) ; [Szilagyi et.al 


(|2000|) ; ICalabrcsc, Lchncr and Tiglioj (p002|) ; |Sziiagyi, Schmidt and Winicom 


Future effort will also be aimed in new directions. For example, it would be 
very desirable to model the late inspiral of binary black holes with a quasi- 
adiabatic approximation, similar to the binary neutron star calculations de¬ 
scribed in Section 10.4. 


Lastly, we have only discussed binaries containing two black holes or two 
neutron stars, while black hole-neutron star binaries have so far been neglected 
(but scc pVlillcrl ( [2001|) ). The inspiral, coalescence and merger of black hole- 
neutron star binaries, including the possible tidal disruption of the neutron 
star, seems like an extremely interesting subject and a promising source of 
gravitational radiation, and is likely to receive more attention in the future. 
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A Notation 


We adopt the notation of [Wald| (|1984|) , which is based on that of [Mistier. 


Thorne and Wheelerj (|1973|). By convention, we also adopt the “fortran” con¬ 


vention in which the letters a — h and o — z are four-dimensional indices and 
run from 0 to 3, while the letters i — n are three-dimensional indices that run 
from 1 to 3. We use geometrized units in which c — G — 1. 


Unless stated differently, objects associated with the four-dimensional met¬ 
ric g a b are denoted with a ^ in front of the symbol, objects associated with 
the conformally related (three-dimensional) metric 7 ^ carry a bar, and only 
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objects related to the spatial metric 7 ^ do not carry any decorations. For 
example, T l jk is associated with 7 ^, f* fc with 77 ,, and with g afe . The co¬ 

variant derivative operator is denoted with D { and Di when compatible with 
the spatial metric and the conformally related metric, but with the nabla sym¬ 
bol V a when compatible with the four-dimensional metric g ab . We occasionally 
use A flat for the flat scalar Laplace operator. 

We denote the symmetric and antisymmetric parts of a tensor with brackets 
() and [], for example 

T(ab) = ~ {T ab + T ba ) and Tj a& ] = ~{T ab — T ba ). (A.l) 


Finally, we write a flat four-dimensional metric as r/ab (Minkowski spacetime) 
and a flat three-dimensional metric as rjij. 


B Solving the Vector Laplacian 


In this appendix we discuss two approaches to solving the flat vectorial Poisson 
equation 

(A l at wy = S\ (B.l) 

or, in Cartesian coordinates, 

d j djWi + ^di&Wj = Si. (B.2) 


We have encountered this operator in the momentum constraint (3.13) and 
(7.2), in the shift condition of the thin sandwich approach (3.28), in the min¬ 
imal distortion condition (5.32), and in the Gamma-freezing condition (5.35). 

Bowen and York ( |1980| ) suggest writing the vector W % as a sum of a vector V t 
and a gradient of a scalar U, 

Wi = Vi + diU. (B.3) 


Inserting this into (B.2) yields 


d j djVi + %d j Vj + d j djdiU + ^<9; d j djU = Si. 

o o 


(B.4) 
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We can now choose U such that the two U terms in (B.4) cancel the second 
Vi term, 


cPd,U = (B.5) 

so that (B.4) reduces to 

PdjVi = Si. (B.6) 


We have hence rewritten the vector Poisson equation (B.2) as a set of four 
scalar Poisson equation for U and the components of Vi. 


A second approach has been suggested by |Shibata| ( |1999c| , see also |Oohara 
Nakamura and Shibata| (|1997|) , Appendix C), who used the ansatz 


w t = V,- l -(d, u + x k 9 t P t ). 


(B.7) 


Inserting this into (B.2) yields 

l&djPi - ^di&djU - ^x k did j djP k = Si (B.8) 

If we now choose U so that it satisfies 

d 3 djU = -S jX j , (B.9) 


then (B.8) reduces to 

l&djPi + ^diSj - ^di&dkPj = |Si (B.10) 

and is solved by 

d j djPi = Si. (B.ll) 


The vector equation (B.2) has again been reduced to a set of four scalar Poisson 
equations, (B.9) and (B.ll). While the approach of [Shibata| ( |1999c|) seems a 
little more complicated, it has the advantage that the source terms in (B.9) and 
(B.ll) are non-zero only where S t is non-zero. In some cases (for example for 
the momentum constraint (3.13)) this may lead to compact sources, which can 
have advantages for numerical implementations (prandclement et. al . |, |2001|). 


110 

























In the approach of [Bowen and York| (|1980|) , on the other hand, the source 
term of equation (B.5) is never compact. A third approach has been suggested 
by IQohara and Nakamura| ( |1999|) . Numerical implementations (using spectral 
methods) of the three approaches have been compared by |Grandclement et. al\ 


C Conformally Flat or Not? 


Many numerical calculations, especially for the construction of initial data, 
assume the spatial metric to be conformally flat. Simultaneously, many au¬ 
thors have pointed out the limitations of conformal flatness, and have argued 
strongly against that simplification. Not all arguments have been correct, and 
it may be useful to briefly review those both in favor and against conformal 
flatness. 


As we have seen in Section 3, conformal flatness greatly simplifies the initial 
value equations. Moreover, before the initial value equations can be solved, 
some form of the conformally related metric has to be chosen. In some cases 
educated guesses can be made (for example by choosing Kerr-Schild back¬ 
ground data ([Marronetti and Matzner[ |2000|) or by adopting a post-Newtonian 
metric instead of a flat conformally-related background metric as suggested in 
Section 7.3). However, in most cases it may not be clear what a better choice 
than conformal flatness might be. 


In Section 3 we found that the dynamical degrees of freedom of the gravi¬ 
tational fields can be identified with parts of the conformally related spatial 
metric and the transverse-traceless part of the extrinsic curvature. This sug¬ 
gests that the assumption of conformal flatness and vanishing of may 
“minimize the gravitational radiation content” of a spatial slice £. This argu¬ 
ment is not strictly true in general; for example it does not even hold for single 
rotating black holes. Rotating Kerr black holes, for example, which are station¬ 
ary and do not emit any gravitational radiation, are not conformally flat 29 . 
Similarly, conformally flat models of rotating black holes that are constructed 
in the Bowen-York formalism do contain gravitational radiation ( [Brandt and 
Scidcl| . |1995a| ,|5. |1996| ; |Gleiser et.al . | , |1998| , see also |Jansen et.al\ ( |2002| )). 


For rapidly rotating single neutron stars it has been shown that conformal 
flatness introduces an error of at most a few percent ( |Cook, Shapiro and 
Tcukolskyl , |1996|) . Similarly small discrepancies were found by [Usui, Uryu ancl 


^ At least slices of constant Boyer-Linquist time are not conformally flat, nor are 
axisymmetric foliations that smoothly reduce to slices of constant Schwarzschild 
time in the Schwarzschild limit flGarat and Price , 2000| ) 
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Eriguchi ( |1999|) ; [Usui and Eriguchi| ( j2002|) , who compared conformally flat 
binary neutron star models with models constructed under different assump¬ 
tions. These small deviations are not surprising, since differences between 
a conformally flat metric and the “correct” metric appear at second post- 
Newtonian order (e.g. |llicth and Schafci] ( |1996|) ), which are the order of a few 
percent for neutron stars. It is therefore quite likely that, at least for neutron 
stars, the error arising from the conformal flatness assumption may be less 
than other errors typically expected in current simulations, including finite 
difference errors and effects from the poor handling of the outer boundaries. 


Wilson and Mathews ( |1995| , see also |Wilsom Mathews and Marronetti| ( |1996|) ) 
adopted the conformal flatness approximation in their simulations of binary 
neutron stars, and found that their neutron stars collapsed to black holes 
individually prior to merging. Since this result was very counter-intuitive, it 
was suspected that this surprising result was erroneous and an artifact caused 
by the assumption of conformal flatness. It was later found that these findings 
were indeed wrong and that they were caused by an error in the derivation 
of the equations ( pianagan[ |1999| ; [Mathews and Wilson[ |2000|) and not by 
conformal flatness. 


It has been argued similarly that the disagreement between numerical ( Cook , 
1994j ; |Baumgartc[ [2000| ) and post-Newtonian (e.g. Parnour, Jaranowski and| 


Schafcij ( |2000|) ) values for the ISCO of binary black holes could be caused by 


the assumption of conformal flatness in the numerical calculations. However, 
the more recent numerical calculations of |Grandclement, Gourgoulhon and 
BonazzoIa| ( [2002| ) achieve much better agreement with the post-Newtonian 
results, even though they also assume conformal flatness (also [Blanclicf| ( |2002| ); 
Damour, Gourgoulhon and Grandcement| ( [2002|) ). As we discussed in Section 
7.3, it is more likely that the choice of initial value decomposition, which 
affects the transverse parts of the extrinsic curvature (see the discussion in 
Pfeiffer, Cook and Tcukolskyl (|2002|) ), caused the earlier discrepancies. The 
good agreement between the numerical results of |Grandclement , Gourgoulhon 


land Bonazzola| ( |2002| ) and post-Newtonian results (|Blanchct| , [2002| ; Parnour 


Gourgoulhon and Grandcement| , |2002| ) may even suggest that the assumption 
of conformal flatness is quite adequate for binary black hole models, at least 
as long as the spin of the individual black holes is not too large. 
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